AD-A118  518  DELAWARE 

numerical 

JUN  82  J 

UNCLASSIFIED  CE-82-24 

UNI  V  NEWARK 
MODELING  OF 

T  KIRBY  *  R  A 

DEPT  OF  CIVIL 
rHE  NEARSHORE 
DALRYMPLE 

ENGINEERING  F/G  8/3 

REGION. (U) 

N00014-81-K-0297 

ML 

fl 

} 

j 

[ 

2 

- - - A 

3 

NUMERICAL  MODELING 
OF  THE  NEARSHORE  REGION 


by 

James  T.  Kir^y,  Jr. 

and 

Robert  A.  Dalrymple 


Technical  Report  No.  11 
Contract  No.  N00014-81-K-0297 

with  the  OFFICE  OF  NAVAL  RESEARCH,  GEOGRAPHY  PROGRAMS 


Research  Report  CE-82-24 
June  1982 

OCEAN  ENGINEERING  PROGRAM 


m 


DEPARTMENT  OF  CIVIL  ENGINEERING 
UNIVERSITY  OF  DELAWARE 
NEWARK,  DELAWARE 
19711 


A 


NUMERICAL  MODELLING  OF  THE  NEARSHORE  REGION 


by 

JAMES  T.  KIRBY,  Jr. 
and 

ROBERT  A.  DALRYMPLE 


Technical  Report  No.  11 
Contract  No.  N00014-81-K-0297 
with  the  Office  of  Naval  Research,  Geography  Programs 


RESEARCH  REPORT  CE-82-24 

June  1982 

Ocean  Engineering  Program 

Department  of  Civil  Engineering 
University  of  Delaware 
Newark,  Delaware 


19711 


security  classification  or  this  pacc  fw**,,  dm*  entered) 


REPORT  DOCUMENTATION  PAGE  befor^cwpleto^forii 


REPORT  NUMBER  2.  GOVT  ACCESSION  HO.  ».  RECIPIENT'S  CATALOG  NUMBER 

ONR  TR.  No.  11 


s.  type  op  report  *  period  covcaco 


4.  TITLE  (end  Subtitle) 


Numerical  Modelling  of  the  Nearshore  Region 


7.  AUTHOR^ 


6.  PERFORMING  ORC.  REPORT  HUM*ER 


B.  CONTRACT  OR  GRANT  NUMBERf*> 


James  T.  Kirby,  Jr.  and  Robert  A.  Dalrymplc  N00014-81-K-0297 


.  PERFORMING  ORGANIZATION  name  and  aodhcss 


»0.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  6  »ORK  UiriT  HUM6ERS 


h.  CONTROLLING  OFFICE  name  and  aooress 


Department  of  Civil  Engineering 
University  of  Delaware,  Newark,  DE  19711 


12.  REPORT  OATS 

June  1982 


13.  NUMBER  OF  PAGES 


I  A.  MONITORING  AGEnCY  name  a  AOORESSf//  different  from  Controlling  Office)  IS.  SECURITY  CLASS,  (of  tMe  report) 


ISA.  DECE  ASSI  Fl  cation/  DOWN  grading 
schedule 


l«.  DISTRIBUTION  STATEMENT  ( o f  i him  Report) 


This  report  has  been  approved  for  public  release  and  sale;  its 
distribution  is  unlimited 


17.  DISTRIBUTION  STATEMENT  (of  the  ab«lr«cl  mntmred  tn  Block  70,  ft  rf///»rml  from  Report) 


IS.  KEY  WORDS  (Continue  on  rereree  mid*  It  n«e memory  end  Identify  by  block  number) 


Numerical  models,  nearshore  circulation,  longshore  currents, 
rip  currents 


20.  ABSTRACT  (Confirm*  on  reeeeee  efde  ((  nm  cm  emery  end  Identify  by  block  number) 

The  thoeretical  background  and  numerical  formulation  for  two 
finite-difference  models  for  predicting  nearshore  circulation  are 
reviewed.  The  models  are  based  on  the  time-  and  depth-averaged 
equations  of  continuity  and  momentum,  with  momentum  flux  terms 
derived  from  linear  wave  theory.  Wave  shoaling  and  refraction  are 
determined  by  the  method  of  Noda  et  al .  (1974).  The  first  model. 


DO  1  j  ah  n  1473  EDITION  OF  f  NOV  «MS  ORSOLEVC 


Unclassified  _ 

SECURITY  CLASSIFICATION  CF  THIS  PAGE  ( 


i  Deim  £nie~e*f) 


*ecu«»ry  classification  o»  th.s  f*ag 


i 


referred  to  as  linear,  neglects  the  effects  of  convective 
acceleration  and  lateral  mixing,  and  calculates  bottom  friction 
using  the  small  mean  current  assumption.  The  second,  nonlinear 
model  includes  the  neglected  effects,  and  uses  the  full  bottom 
stress  relation  as  given  by  Liu  and  Dalrymple  (1978) .  Both 
models  solve  iteratively  for  mean  currents  and  mean  free  surface 
displacement  at  discrete  grid  points. 

A  calibration  of  both  models  is  described  based  on  comparison 
with  field  data  obtained  from  the  NSTS  Torrey  Pines  experiment 
(Gable,  1979) .  Finally,  various  results  obtained  during  previous 
studies  utilizing  the  models  are  presented  and  discussed. 


Unclassified _ 

SECURITY  CLASSIFICATION  OF  THIS  ^AGEf"***  *>«*• 

ii 


ABSTRACT 


The  theoretical  background  and  numerical  formulation  for  two 
finite-difference  models  for  predicting  nearshore  circulation  are  reviewed. 
The  models  are  based  on  the  time-  and  depth-averaged  equations  of  continuity 
and  momentum,  with  momentum  flux  terms  derived  from  linear  wave  theory.  Wave 
shoaling  and  refraction  are  determined  by  the  method  of  Noda  et_  al^.  (1974). 
The  first  model,  referred  to  as  linear,  neglects  the  effects  of  convective 
acceleration  and  lateral  mixing,  and  calculates  bottom  friction  using  the 
small  mean  current  assumption.  The  second,  nonlinear  model  includes  the 
neglected  effects,  and  uses  the  full  bottom  stress  relation  as  given  by 
Liu  and  Dalrymple  (1978).  Both  models  solve  iteratively  for  mean  currents 
and  mean  free  surface  displacement  at  discrete  grid  points. 

A  calibration  of  both  models  is  described  based  on  comparison  with 
field  data  obtained  from  the  NSTS  Torrey  Pines  experiment  (Gable,  1979). 
Finally,  various  results  obtained  during  previous  studies  utilizing  the 
models  are  presented  and  discussed. 


i 


TABLE  OF  CONTENTS 


Page 

Abstract  .  i 

List  of  Figures .  vii 

Chapter  I.  Introduction . . .  1 

Chapter  II.  Theoretical  Development  of  the  Models  .  3 

1.  Introduction  .  3 

2.  Theoretical  Framework  of  the  Models  .  3 

2.1  Specification  of  the  Boundary  Value  Problem  ...  5 

2.2  Depth  and  Time  Averaged  Forms  of  the  Equations  .  .  9 

2.3  Radiation  Stresses  .  12 

2.4  Wind  Stress  . .  13 

2.5  Bottom  Stress .  14 

2.6  Wave  Refraction  and  Shoaling  Including 

Wave-Current  Interaction  .  15 

2.7  Wave  Breaking  Criteria  .....  .  19 

2.8  Lateral  Mixing . 20 

2.9  Wave  Height  Decay .  22 

3.  Overview  of  the  Linear  and  Nonlinear  Models . 25 

Chapter  III.  Formulation  of  the  Numerical  Models  .  29 

1,  Introduction .  29 

2,  The  Grid  Scheme  and  Location  of  the  Unknown 

Quantities .  30 

3,  Boundary  Conditions  .  34 


Page 

4.  Finite  Difference  Operators  and  Notation  .  36 

5.  Finite  Difference  Forms  of  the  Governing 

Equations  -  Linear  Model  .  37 

6.  Finite  Difference  Forms  of  the  Governing 

Equations  -  Nonlinear  Model  .  40 

7.  The  Numerical  Scheme  for  Refraction  and  the  Wave 

Height  Field  .  43 

Chapter  IV.  Calibration  of  the  Nearshore  Circulation  Models  ...  49 

1.  Field  Data  Used  for  Calibration . 49 

1.1  Choice  of  Field  Data  Used  for  Calibration  ....  51 

2.  Calibration  of  the  Models . 61 

2.1  Linear  Model  Calibration  .  63 

2.2  Nonlinear  Model  Calibration  .  65 

2.3  Response  of  Both  Models  Using  Data  Set  2 . 68 

3.  Discussion  of  the  Calibrated  Coefficients  .  72 

Chapter  V.  Examples  of  Numerical  Results  Using  the  Nearshore 

Circulation  Models  .  77 

1.  Introduction .  77 

2.  Specific  Applications  of  the  Linear  and  Nonlinear 

Models .  77 

2.1  Plane  Beach  Applications  .  78 

2.2  Barred  Profile  Applications  .  86 

2.3  Applications  to  the  Laboratory  Wave  Basin  ....  90 

2.4  Periodic  Bottom  Topography  .  91 

2.5  Intersecting  Waves  Applications  .  102 

Chapter  VI.  Conclusions  .  Ill 


iv 


References 


Appendix  I.  Using  the  Nearshore  Circulation  Program 

Listing  of  Linear  Model  . 

Listing  of  Nonlinear  Model  . 


v 


LIST  OF  FIGURES 


Figure 

No .  Title  Page 

2-1  Plan  Definition  Sketch  For  Coordinate  System  6 

2- 2  Longue t-Higgins f  Analytical  Solution  for  Oblique 

Wave  Attack  on  a  Plane  Beach  21 

3- 1  Grid  Scheme  31 

3-2  Velocity  Components  for  Grid  Block  (i,j)  33 

3- 3  General  Leapfrog  Solution  Scheme  43 

4- 1  Location  of  NSTS  Torrey  Pines  Experiment  50 

4-2  Instrumentation  for  NSTS  Torrey  Pines  Experiment  52 

4-3  Torrey  Pines  Nearshore  Bathymetry,  November  9,  1978  53 

4-4  Torrey  Pines  Nearshore  Bathymetry,  November  18,  1978  54 

4-5  Time  Averaged  Current  Measurements,  November  10,  1978  57 

4-6  Average  Field  Velocity  Profile,  Data  Set  1  58 

4-7  Angular  Distribution  of  Spectral  Density,  November  15, 

1978  "  60 

4-8  Average  Field  Velocity  Profile,  Data  Set  2  62 

4-9  Currents  Induced  in  Model  Using  Nov.  9  Bathymetry. 

Nonlinear  Model :  Data  Set  1  64 

4-10  Longshore  Average  Depth  Profile,  Torrey  Pines,  Nov.  9  66 

4-11  Variation  of  Longshore  Current  with  Friction  Factor 

in  the  Linear  Model  67 

4-12  Variation  of  Longshore  Current  with  Lateral  Mixing  in 

the  Nonlinear  Model.  f  =  0.015  69 

4-13  Variation  of  Longshore  Current  with  Lateral  Mixing  in 

the  Nonlinear  Model.  f  =  0.010  70 

4-14  Longshore  Current  in  the  Linear  and  Nonlinear  Models. 

Data  Set  2  71 


vii 


K  igure 

No .  Ti  1 1 1.-  P.age 

-4-15  Currents  Induced  in  Model  Using  Nov.  18  Bathymetry. 

Nonlinear  Model:  Data  Set  2  73 

5-1  Set-up  on  a  Plane  Beach  79 

5-2  Set-up  at  Shoreline  Duo  to  a  Wave  Croup  with  a  Period 

of  18  Seconds  81 

5-3  Resonance  of  Wave  Channel  Du«.  to  Forcing  at  Seiche 

Period  82 

5-4  Inshore  Mean  Free  Surface  Displacement  Versus  Time 

for  the  Nonlinear  Model  Application  to  a  Plane  Beach  83 

5-5  Longshore  Current  Profile  for  Plane  Beach. 

Linear  Mode  1  8 4 

5-6  Longshore  Current  Profile  for  Plane  Beach. 

Nonlinear  Model  85 

5-7  Btirred  and  Plane  Beach  Profiles  87 

5-8  Longshore  Current  Profile  for  Barred  Beach.  Linear 

Model  88 

5-9  Longshore  Current  Profile  for  Barred  Beach.  Nonlinear 

Model  89 

5-10  Experimental  Streamlines  in  the  Wave  Basin  (From 

Dalrymple  e_t  <^1.  ,  1977)  92 

5-11  Currents  in  Model  Basin  Corresponding  to  Dalrymple  et_  al . 

(1977)  '  93 

5-12  Predicted  and  Measured  Longshore  Velocities  for  the  Basin 

Test  Case  94 

5-13  Depth  Contours  for  the  Periodic  Bottom  of  Noda  (1973)  96 

5-14  Current  Vectors  for  Noda  Topography.  Linear  Model  97 

5-15  Current  Vectors  for  Noda  Topography.  Nonlinear  Model  98 

5-16  Currents  Induced  by  Bulge  on  a  Plane  Beach.  Linear 

Model  100 

5-17  Currents  Induced  by  Bulge  on  a  Plane  Beach.  Nonlinear 

Model  101 


vi  i  i 


F igure 
No. 


T  i  1 1  e 


.p-aj.c‘ 


1 

cc 

Definition  Sketch  for  the  Intersecting  Wave 

App 1 icat ion 

102 

5-19 

Current  Vector  Plot  for  a  Rip  Current  Perpendicular 
to  the  Shoreline 

108 

5-20 

Current  Vector  Plot  tor  the  Meandering  Circulation 
Pattern 

109 

‘ 


Chapter  I 


INTRODUCTION 


The*  need  for  methods  to  accurately  predict  the  magnitude  and  spatial 
distribution  of  nearshore  currents  is  central  to  the  present  research  efforts 
aimed  at  quantifying  the  transport  of  marine  sediment  in  the  beach  and  near¬ 
shore  environment.  Investigators  have  made  significant  strides  in  describing 
the  mean  wave- induced  motions  in  the  surf  zone  in  the  twenty  years  since  the 
correct  formulation  of  the  averaged  equations  of  motion.  Using  the  concept 
of  wave  momentum  flux,  or  radiation  stress  (see  Longuet-Higgins  and  Stewart, 
1964),  it  has  become  possible  to  analytically  predict  the  wave  set-up,  long¬ 
shore  currents,  and  spacing  of  rip  currents  on  beaches  of  simple  planfom. 

Recent  strides  have  also  been  made  towards  predicting  the  dynamic  response 
of  the  surf  zone  to  fluctuating  driving  forces  with  time  scales  larger  than 
that  of  the  incident  wind  waves,  such  as  surf  beat  (see,  for  example,  Symonds, 
Huntley,  and  Bowen ,  1982). 

In  general,  however,  the  nearshore  environment  is  a  complex  system 
which  is  not  amenable  to  analytic  treatment.  Even  in  simple  situations, 
the  presence  of  a  physical  feature  predicted  by  one  mechanism,  such  as  the 
longshore  current,  will  greatly  complicate  the  prediction  of  a  separate  feature, 
such  as  rip  currents.  In  addition,  the  action  of  waves  and  currents  on  the 
bottom  can  create  complex  topographies,  making  a  reduction  of  the  analytic 
problems  to  one  space  dimension  impossible. 

This  report  represents  a  review  and  conclusion  of  several  studies 
conducted  at  the  University  of  Delaware  with  the  aim  of  providing  a  numerical 
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mo do 1  for  calculating  nearshore  wave-induced  currents  and  mean  water  level 
l luc t ua t ions .  The  purpose  of  constructing  a  numerical  model  rests  on  the 
need  to  extend  our  predictive  capabilities  into  situations  which  lie  beyond 
the  scope  of  analytic  methods.  In  the  end,  all  numerical  models,  as  well 
as  analytic  formulations,  are  limited  in  scope  by  the  simplifying  assumptions 
incorporated  in  their  theoretical  framework;  in  this  regard,  the  present  models 
represent  an  attempt  to  extend  present  analytic  treatments  to  the  case  of  a 
complex  topography  in  two  dimensions.  The  models  do  not  consider  the  associated 
sediment  transport  problem,  although  this  capability  can  be  added  (see  McDougal, 
1979,  and  Paddock  and  Ditmars,  1981).  Also,  the  models  require  that  the 
incident  wave  field  be  regarded  as  monochroma t ic ,  or,  after  some  model 
modifications,  narrow  banded  enough  to  be  represented  as  a  modulated  wave 
train  at  a  single  carrier  frequency. 

The  models  described  here  utilize  a  wave  refraction  scheme  developed 
by  Noda  e_t  a_l.  (1974).  Using  this  scheme,  Birkemeier  and  Dalrymple  (1976) 
developed  a  circulation  model,  referred  to  here  as  the  linear  model,  which 
neglected  the  effect  of  convective  accelerations  and  lateral  mixing.  The 
model  of  Ebersole  and  Dalrymple  (1979),  referred  to  here  as  the  nonlinear 
model,  extended  the  treatment  to  include  these  effects. 

In  Chapter  II,  the  theoretical  framework  for  the  two  models  is 
described,  followed  by  an  outline  of  the  numerical  formulations  in  Chapter 
III.  In  Chapter  IV,  we  present  a  calibration  of  the  models  using  field  data 
from  the  NSTS  Torrev  Pines  experiment  (Gable,  1979).  Chapter  V  gives  a 
summary  of  appl icat ions  of  the  models  presented  in  Birkemeier  and  Dalrymple 
(1976)  and  Ebersole  and  Dalrymple  (1979). 
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Chapter  II 


THEORETICAL  DEVELOPMENT  OF  THE  MODELS 


1.  INTRODUCTION 

In  this  chapter  the  theoretical  framework  for  the  development  of 
time  averaged  governing  equations  for  the  problem  of  waves  and  currents  in 
the  nearshore  zone  is  outlined.  The  development  of  the  numerical  circu¬ 
lation  models  in  either  the  ’’linear"  or  "nonlinear"  form  is  then  described 
in  Chapter  III  based  on  the  theoretical  framework. 

In  general,  the  development  of  each  model  has  been  described  in 
previous  technical  reports;  the  linear  model  in  Birkemeier  and  Dalrvmple 
(1976),  and  the  nonlinear  model  in  Ebersole  and  Dalrvmple  (1979).  For  this 
reason,  some  of  the  derivations  and  intermediate  steps  needed  to  develop 
the  governing  equations  are  not  described  in  detail  in  the  present  report. 

The  reader  can  refer  to  the  previous  work  for  missing  details.  However, 
both  models  now  include  the  option  of  calculating  wave  energy  decay  due  to 
interaction  with  the  bottom,  which  has  not  been  included  in  previous  reports. 
The  theory  and  implementation  of  this  option  is  described  in  detail. 

2.  THEORETICAL  FRAMEWORK  FOR  THE  MODELS 

The  basis  for  any  fluid  dynamic  model  rests  on  the  principle  of 
conservation  of  mass,  conservation  of  momentum  (the  Navier-Stokes  equations), 
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and  conservation  of  energy.  The  resulting  system  of  equations,  together 
with  boundary  conditions  which  quantify  the  interaction  o!  the  fluid 
continuum  with  its  solid  and  free  bounding  surfaces,  give  a  mathematical 
representation  of  the  physical  problem  of  interest.  The  problem  may  then 
be  further  simplified  by  assumptions  which  are  consistent  with  the  physical 
processes  involved. 

Here,  we  are  interested  in  the  effect  of  waves  propagating  towards 
shore  over  a  complex  bathymetry  and  breaking,  and  the  mean  currents  driven 
by  changes  in  the  flux  of  wave  momentum.  The  problem  has  two  apparent  time- 
scales;  a  fast  time  scale  associated  with  the  oscillation  of  the  incident 
waves,  and  a  slower  time  scale  over  which  the  characteristics  ol  the  incident 
wave,  such  as  height  and  angle  of  incidence,  may  vary.  The  longer  time  scale 
may  also  include  the  effect  of  changes  in  wind,  tidal  oscillation,  and,  in 
models,  which  include  sediment  transport,  gradual  shift  of  the  bottom.  Since 
our  attention  here  is  towards  mean  quantities  which  are  reasonably  steady 
in  time,  the  set  of  equations  may  be  averaged  over  the  faster  time  scale  of 
the  wave  oscillation  to  remove  the  direct  effect  of  the  oscillation.  Tin- 
effect  of  the  waves  then  is  represented  by  a  stress- like  term  acting  on  the 
slowly  varying  mean  flow  pattern.  In  addition,  since  we  are  mainly  interested 
in  net  transport  quantities  rather  than  detailed  structure  of  the  velocity 
profiles  over  depth,  the  equations  may  be  averaged  over  depth,  reducing  the 
entire  problem  to  a  two  dimensional  problem  in  the  horizontal  plane  together 
with  appropriate  boundary  conditions.  This  averaged  model  can  then  be  solved 
using  numerical  procedures,  as  discussed  in  Chapter  III.  The  quantities  to 
be  determined  are:  the  horizontal  components  of  mean  wave-induced  current. 
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the  local  wave  height  and  angle  of  incidence,  and  the  set-up,  or  wave-induced 
deviation  of  the  mean  water  surface  from  its  still-water  level. 


2 . 1  Specification  of  the  Boundary  Value  Problem 


A  right-handed  coordinate  system  is  defined  with  x  in  the  offshore 
direction,  normal  to  the  coastline,  y  in  the  longshore  direction,  and  z 
vertically  upward  (see  Figure  1).  The  continuity  equation  is 

3p  .  3(pu)  3 (pv)  3(pw)  _  / 

at  ax  ay  az  v 

where  p  is  the  water  density  and  (u,v,w)  are  the  (x,y,z)  components  of 
velocity,  respectively.  The  continuity  equation  can  be  further  reduced 
to  the  form 


3u  3v  3w 
Dx  +  +  az 


0 


(2.2) 


consistent  with  an  assumption  of  constant  water  density  p. 


The  momentum  equations  are,  in  the  x  direction 

_  „  2  „  .  „  ,  3t  3x  3t 

9u  3u  9uv  3uw  __  ^  1  DP  ^  1  ,  xx  ,  yx  zx . 

at  3x  3y  9z  p  3x  p  3x  3y  3z 

in  the  v  direction, 


_  _  ^2  ^  n  3x  3t  3t 

av  3uv  3y  3vw  _  1  3P  ,  1  r  xy  yy  +  zy  ■> 

3 1  3x  9y  3z  p  3y  p  3x  3y  3z 


(2.4) 


and  in  the  z  direction, 


3w  9uw  9vw  3w‘ 
3t  3x  3y  3z 


1  _3_ 
p  3z 


(P  +  gz) 


3t  9t  3t 
xz  yz  +  zz 

3x  3y  3z  ’ 


(2.5) 
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at  ter  substitution  of  the  continuity  equation  into  the  convective  acceleration 
terms . 


WIND  ANGLE 


WAVE  ANGLE 


Figure  2-1. 


Plan  Definition  Sketch  For  Coordinate  System. 


Boundary  Conditions 


Certain  boundary  conditions  are  required  at  the  physical  boundary 
of  the  water  body  in  question  in  order  to  correctly  specify  the  problem. 
First,  kinematic  conditions  are  specified  at  the  free  surface  and  rigid 
bottom,  which  state  that  water  particles  may  not  cross  the  boundary 
surface,  whether  it  be  rigid  or  moving. 

The  equation  of  a  surface  of  the  fluid  domain  is  given  by 

F(x,y,z,t)  =  0 

A  water  particle  cannot  flow  across  the  surface,  otherwise  the  surface 
wouLd  cease  to  exist.  Mathematically,  this  condition  is  expressed  by  the 
total  time  rate  of  change  of  the  function  F(x,y,z,t)  being  equal  to  zero. 

(F(x,y,z,t))  =  0 

At  the  free  surface  the  boundary  is  given  by 
F^ (x,y , z  ,  t)  =  z  -  n(x,y,t)  =  0 

and  at  the  bottom 
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For  the  bottom  boundary  condition  (BBC)  we  get. 


3h 

3t 


+ 


-h 


_3h 

3x 


+  v 


-h 


dh 

3y 


+  W 


0 


(2.7) 


where  u,v,w  are  the  velocity  components  in  the  x,y  and  z  directions  and  the 
subscripts  denote  the  evaluation  of  a  specific  term  at  the  bottom,  z  =  -h, 
or  at  the  free  surface,  z  =  n . 


In  addition  to  the  kinematic  free  surface  boundary  condition,  a  dynamic 
free  surface  boundary  condition  is  required  as  well,  which  states  that 

P  =  constant  z  =  n 


This  condition  is  satisfactory  if  we  neglect  the  local  generation 
of  waves  by  wind  or  the  deformation  of  the  free  surface  due  to  barotropic 
effects.  Then,  using  Bernoulli’s  equation  expanded  to  the  free  surface, 
we  obtain 


?  2  2 
3(f>  (u"  +  v  +  w  ) 

n  ,  n  n  n 

3 1  z 


+  gn  =  0 


z  =  n 


(2.8) 


after  setting  P  =  0,  where  <)>(x,y,z,t)  is  the  velocity  potential  for  the  wave 
motion . 


The  model  also  requires  lateral  boundary  conditions.  In  the  y 
(longshore)  direction,  the  bathymetry,  given  by  the  surface  h(x,v)  will  be 
required  to  be  periodic.  Since  the  offshore  wave  conditions  will  be  assumed 
uniform,  lateral  periodicity  conditions  for  wave  and  currents  are  also  assumed. 

Offshore  in  the  x  direction,  the  usual  boundary  condition  for  a  wave 
problem  would  be  to  assume  that  all  waves  at  the  boundary  other  than  the 
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incident  wave  are  propagating  away  from  shore ,  the  radiation  condition 
(Sommerfeld  (1949)).  However,  the  circulation  model  does  not  directly 
calculate  an  actual  wave  field.  Rather,  a  somewhat  arbitrary  condition, 

u  =  0  ;  x  =  furthest  offshore  grid 

is  imposed,  which  serves  to  put  a  bound  on  the  horizontal  extent  of  the  flow 
under  consideration. 

Specification  of  an  onshore  boundary  condition  is  an  uncertain  task, 
due  to  the  complexity  of  the  surf  zone.  In  general,  it  is  likely  that  some 
wave  energy  survives  the  breaking  process  and  reflects  from  the  beach,  leading 
to  waves  propagated  back  into  the  region  of  the  model.  However,  it  is  assumed 
for  the  purpose  of  modelling  that  all  wave  energy  decays  in  the  surf  zone, 
reducing  the  wave  height  to  a  value  of  zero  at  the  shoreline.  In  addition , 
both  u  and  v  are  set  equal  to  zero  at  the  shoreline. 

2 . 2  Depth  and  Time  Averaged  Forms  of  the  Equations 

By  integrating  the  equations  of  motion  and  continuity  over  depth  and 
substituting  the  boundary  conditions,  the  problem  is  reduced  to  equations  in 
two  horizontal  dimensions,  with  lateraJ  boundary  conditions  prescribed. 
Secondly,  the  quantities  of  principal  interest  are  average  in  nature,  i.e., 
mean  currents,  wave  height  and  wave  angle,  and  mean  water  level.  The  equations 
can  then  be  time  averaged  over  a  wave  period  to  remove  consideration  of 
instantaneous  wave  induced  motions.  The  quantities  remaining  would  be  free 
to  vary  slowly  in  time  in  response  to  changing  offshore  conditions. 
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Integrating  Eq .  (2.2)  over  depth  from  z  =  -h(x,y,t)  to  z  =  q(x,y,t), 
using  Leibnitz  rule  (Hildebrand  (1976),  p.  365)  to  remove  partial  derivatives 
from  within  the  integrals,  and  substituting  Eqs.  (2.5)  and  (2.6),  the 
continuity  equation  becomes 


_3_ 

3t 


(2.9) 


Let  both  u  and  v  be  comprised  of  a  time  independent  mean  flow  and 
a  wave  induced  flow, 

u  =  U  +  u 
v  =  V  +  v 


By  substituting  the  above  expressions  for  u  and  v  into  Eq .  (2.7) 
and  time  averaging  over  one  wave  period  so  as  to  eliminate  the  wave  induced 
fluctuations,  the  continuity  equation  can  be  written  as 


it {p(h  +  n)}  +~h  {pU(h  +  n)}  +~k 


a  —  —  3  r 

+  —  {pV (h  +  q)  }  +  —  pvdz  =  0  (2.10) 

dy  J_h 

where  the  symbol  " -  "  denotes  the  time  average  over  one  wave  period  and 

n  the  time  independent  mean  free  surface  displacement.  Note  that  the  time 
averages  of  the  vertically  integrated  wave  induced  velocities  u  and  v  are  not 
zero . 

Defining  U  =  U  +  u 

V  H  V  +  v 
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pud 2  pvdz 

t  -  J-h  ,  J-h 

where  u  =  - ~  and  v  =  - - — 

p(h  +  n)  p  (h  +  n) 

are  the  wave  induced  mass  transport  velocities,  and  substituting  the  total 
depth  D  for  (h  +  n) ,  the  time  averaged,  depth  integrated  continuity  equation 
is,  in  its  final  form 

It  -  h  <UD)  +  i  <VD>  ■ 0  •  (2-n> 

Inherent  in  the  derivation  are  the  assumptions  that  the  bottom  is  constant 
with  time  and  the  density  is  constant  in  space  and  time. 


Momentum  Equations 


The  x  and  v  momentum  equations  are  manipulated  in  the  same  way  as  the 
continuity  equation  in  order  to  achieve  equations  which  are  independent  of  wave 
induced  oscillations,  i.e.,  they  are  integrated  over  depth  and  time  averaged 
over  a  wave  period.  Details  of  this  derivation  may  be  found  in  Ebersole  and 
Dalrvmple  (1979) .  The  resulting  equations  are  equivalent  to  those  given  by 
Phillips  (1977),  but  are  written  explicitly  in  terms  of  depth  averaged  mean 
velocities.  The  equations  contain  terms  for  mean  bottom  shear  stress,  mean 
surface  shear  stress  due  to  wind,  mean  lateral  friction  (not  used  in  the 
linear  model),  and  excess  mean  momentum  stress  due  to  wave  action  (the 
radiation  stress,  see  Longuet-Higgins  and  Stewart  (1964)).  The  x  momentum 
equation  in  its  final  form  can  be  written  as, 


it  (UD)  + 


8x 


(u2d) 


d_ 

dv 


(UVD)  =  -gD 


2p 

3x 


3v 


,  9S  1  3S  .  , 

1  XX  1  XV  .  1 - 1  - 

—  - - - — +  —  T  “  —  T, 

p  3x  P  3y  p  sx  p  bx 


(2.12) 


and  the  final  form  of  the  y  momentum  equation  can  he  written  as 


3t|  D  ■>TC. 


57  (VD)  +  (UVD)  +  f  (V2D)  =  -gD 

ax  Jv  3v  o  d 


p  3x 


i  i  3S  i  i 

-I  xy  _  1  y  v  1 - JL_  - 

o  3x  o  3v  p  sv  p  bv 


(2.13) 


Wave  Energy  Equation 

Following  Phillips  (1977),  an  equation  expressing  the  conservation  of 
averaged  wave  energy  may  he  written  as 

+  d-  {E(U  +  C  cos  0) i  +  ~  {E(V  +  C  sin  0)} 

It  lx  g  3y  g 


„  3U  (3V  3V 

+  s  —  +  S  r- - F  +  S  =  e 

xx  Ox  xv  Ox  Ov  vv  Ov 


(2.14) 


Here,  C  is  the  wave  group  velocity  and  4  is  the  local  wave  angle.  The 
quantity  e  represents  the  dissipation  of  wave  energy,  and  is  identically  zero 
in  a  conservative  wave  field.  It  can  be  given  a  non-zero  value  to  include  the 
effect  of  wave  damping,  as  described  below. 


2 . 3  Radiation  Stresses 

The  radiation  stresses  included  in  Eqs.  (2.12),  (2.13)  and  (2.14) 
represent  the  stress  on  the  water  column  induced  by  wave  action.  Neglecting 
the  effects  of  small  scale  turbulent  velocities,  the  stresses  have  been  given 
in  simple  form  by  Longuet-H xggins  and  Stewart  (1964)  as 

S  =  E  [(2n  -  l/2)cos^0  +  (n  -  l/2)sin"01  (2,15) 


Svv  -  E  I(2n  -  l/2)sin“6  4-  (n  -  l/2)cos^6]  (2.16) 

SXV  =  svx  =  !  n  8ln  (29)  (2-1 7) 


where 

E  is  tht 

*  wave  ene 

rgy 

(C  )  , 

to  wave 

celerity 

(C) 

wave , 

E  and  n 

are  given 

by 

6  is  the  wave  angle,  and  n  =  ratio  of  group  velocity 
To  second  order  for  a  progressive  small  amplitude 


E 

n 

k 

h 

L 

H 


"  8  P8  H‘ 


C 


-  1.  [I  4-  - 2kh - j 

C  2  sinh(2kh)‘ 


=  wave  number  (= 
~  depth 
=  wave  length 
=  wave  height. 


(2.18) 

(2.19) 


2 .4  Wind  Stress 

Although  other  methods  exist  for  computing  the  surface  stress  due  to  the 
wind  (see  Wu  (1968)),  the  one  suggested  in  the  Shore  Protection  Manual  (1977) 
has  been  utilized.  This  form  was  first  developed  by  Van  Dorn  (1953)  and  gives 
a  fairly  good  fit  to  the  existing  experimental  data.  The  form  of  the  surface 
stress  is  quadratic  in  the  wind  speed  and  is  given  by 

t  =  pKlwjw  (2.20) 

sx  '  '  X 

I 

t  =  ptdwlw  (2.21) 

sy  v 

where  W  is  the  wind  speed,  and  W  ,  arc  wind  velocity  components  in  the  x  and 
y  directions,  as  determined  hy  the  wind  angle,  a. 

I 
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TIk-  wind  stress  coe  f  i  i  e  lent  K  is  determined  em;>  i  r  i  ca  1  1  v 
dependent  on  the  magnitude  of  the  wind  Velocity  sueh  that 

K, 


to  he 


and 


K  - 


/ 

W 


1 


K.  +  K  ,  (1  +  K  /W)  ‘ 
1  2  cr 


K  j  =  1.1  x  10 


-b 


2.5  x  10 


-6 


(2.22) 


W  =  7.2  mu ter s/ second 

c  r 

A  comparison  of  the  wind  stress  coefficient  given  here  with  experimental 
data  is  given  in  Pearce  (1972). 

2.5  Bottom  Stress 


The  correct  method  for  specifying  the  effect  of  a  rigid  bottom  on 
waves  and  currents  is  still  a  matter  of  lively  debate.  For  the  purpose  of 
modelling  hydrodynamics  in  the  nearshore  zone,  the  average  bottom  shear  stress, 


,  is  generally  taken  to  be  of  the  form 


Tb  =  p  8  vv 


(2.23) 


where  f  is  a  Darcv-Wei sbach  friction  factor  and  u,  is  the  total  instantaneous 

b 

scalar  velocity  at  the  bottom  (Longuet-Higgins  (1970a)).  This  relation  is 
known  as  a  quadratic  friction  law.  The  components  of  shear  stress  in  the  x 
and  y  directions  can  be  given  quite  generally  by  (Liu  and  Dal rymple  (1978)) 


_  f  fZ7T  _js 

t,  =  77 —  (U  +  u  cos  0  cos  at)  *  j  u,  |  d  (at) 

bx  16tt  Jq  m  b 


t t  -  77-  I  (V  +  u  sin  0  cos  at)  •  (u.  |  d (at) 
by  16tt  J  q  m  b 


2tt 


(2.24) 


(2.25) 


1  4 


whore  u  is  the  maximum  wave  orbital  velocity  given  i>v 
m  ' 


H 


km  2  sinh  kh 


(2.26) 


and  u ,  is  the  vector  velocity  at  tilt*  bottom.  The  nonlinear  model  of  F.bersoK 
b 

and  Dalrymple  (1979)  retains  these  forms,  whore  the  integration  is  approxi¬ 
mated  by  a  16  term  Simpson's  rule  summation. 

The  linear  model  retains  the  form  originally  used  by  Birkemeier  and 
Dalrymple  (1976).  For  this  form,  the  assumption  that  friction  is  primarily 
due  to  tiie  influence  of  the  wave  orbital  velocity  is  used.  Making  a  small 
mean-current  assumption,  LeBlond  and  Tang  (1974)  show  that  the  bottom  stress 
is  anisotropic,  with 


- _  Pf 

T ,  r  =  —  U  l; 

bx  2:t  m 

T~  =  ftu  v 

bv  4~r  m 


(2.27) 

(2.28) 


with  U,  V,  u  ,  and  f  as  previously  defined.  A  derivation  of  these  equations 
m 

is  given  in  Birkemeier  and  Dalrymple  (1976),  Appendix  A. 


2 . 6  Wave  Refraction  and  Shoaling  Including  Wave-Current  Interation 

The  equations  which  govern  both  wave  refraction  and  shoaling  as  a 
result  of  wave-current  interaction  used  in  the  model  are  those  of  Noda  et  al  . 
(1974).  The  advantage  of  Noda T s  method  is  that  it  can  predict  the  wave  angles 
and  wave  heights  at  certain  points  rather  than  along  a  wave  ray.  This  procedure 
lends  itself  well  to  use  in  the  finite  difference  model  because  calculations 
are  performed  at  points  which  lie  in  the  center  of  rectangular  grid  elements 
which  are  part  of  a  larger  grid  mesh. 


r 


Starting  with  a  progressive  linear  gravity  wave,  the  free  surface 
can  be  written  as, 

■  ;(x,y,t)  =  a (x, y , t ) cos{ 9  (x,v ,  t )  } 

where  a  is  the  wave  amplitude  and  y  is  a  phase  function.  A  wave  number  vector 
can  he  defined  as 

k  -  7*  (2.29) 


and  a  wave  scalar  frequency  can  be  defined  as 


(2.30) 


Using  the  mathematical  property  that  the  curl  of  a  gradient  is  identically 
zero,  it  is  show  that 

7  x  V  =  0 


which  implies  that 

7  x  k  =  0 


This  equality  states  that  the  wave  number  vector  is  irrotat ional .  Assuming 
4>(x,y,t)  is  continuous,  then 

• 

On  substituting  Eqs.  (2.29)  and  (2.30)  into  the  above  expression,  it  is  found 
that 

+  Vo  =  0  (2.31) 

which  is  the  classical  conservation  of  waves  equation. 
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For  tli 


cast*  of  a  wave  propagating  on  a  current  with  velocity 


u  -  ui  +  vj  ,  it  can  bo  shown  that  the  scalar  frequency  with  respect  to 
a  stationary  rtf  crones  franc-  is 


s 

=  a  +  k  •  U 


The  wave  frequency  with  respect  to  a  moving  reference  frame  is  given  by 
the  dispersion  relation, 

a*~  ~  gk  tanh  kh  .  (2.32) 


If  it  is  also  assumed  that  the  wave  number  field  is  constant  in  time  then, 
f rom  Eq .  (2.31), 

V (o  +  k  •  U)  =  0 


or 

- y 

o  +  k  •  U  =  constant 


(2.33) 


This  constant  can  be  evaluated  for  the  case  where  U  =  0  in  which  case 
2tt 

o  =  ~  where  T  is  the  wave  period.  Eq.  (2.33)  then  becomes 


0  4*  k 


(2.34) 


Using  the  coordinate  system  shown  in  Figure  1  and  expanding  Eqs.  (2. 
and  (2.34)  into  Cartesian  coordinates  and  using  the  dispersion  relation,  the 
equations  which  govern  wave  refraction  through  wave-current  interaction  are 
given  by 


r 30  1  Bkx  ,  ,  Q  fB6  ,  1  3k,  n 

cos  0  {- - 7“  — }  4-  sin  0  {—  4-  —  — }  -  0 

3x  k  3y  3y  k  3x 


(2.35) 


29) 
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(2.36) 


igk  tanh(kh)J^ta+  Uk  cos  0  +  Vk  sin  0  =  r.y 


where  0,  k,  h,  U  and  V  are  ail  functions  that  may  vary  in  both  the  x  and  v 
direct  ions. 


The  shoaling  of  waves  is  also  affected  by  the  interaction  of  waves 
and  currents.  The  effect  on  the  waves  is  determined  by  solving  the  energy 
equation.  Dividing  Eq .  (2.14)  by  E  and  expanding  in  Cartesian  coordinates, 
we  get 


—  +  (U  +  C  cos  6)  i  +  (V  +  C  sin  0)  -p  ~ 
E  3t  g  E  3x  g  E  dv 


+  4-  <U  +  C  cos  ©)  +  ~  (V  +  C  sin  0) 

3x  g  3v  g 

+  1{S  f+s  ^  +  s  +  s  2Y>-f- 

E  xx  3x  xy  3y  yy  3y  xy  3v  E 


Using  this  result,  carrying  out  the  differentiation,  and  letting  a 
quantity  Q  be  defined  as 

n  -  i  rq  9U  IV  3V, 

^  E  xx  3x  xy  3y  xy  3x  yy  3x 


the  energy  equation  becomes, 

coS  e)|f  +  (v  +  cgSin  +  f  +  f 

3C  3  c 

-  C  sin  0  —  +  cos  0  -r — ^  +  C  cos  0  +  sin  0  +  Q  -  fr  .  (2.37) 

g  3x  3x  g  3y  3y  E 

For  all  applications  of  the  model  the  wave  height  H  is  assumed  constant  in 
3H 

time,  so  —  =  0.  From  linear  wave  theory  the  group  velocity  C  is  given  by 
3 1  & 
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where 


C 

8 


c 


+  _ } 

sinh(2kli) 


C  -  tanh(kh)  } 

is  the  wave  speed  or  celerity,  k  is  the  wave  number,  and  h  is  the  water  depth. 

2 . 7  Wave  Breaking  Criteria 

Since  Eq.  (2.37)  is  applicable  only  in  determining  the  wave  heights 
of  nonbreaking  waves,  some  method  is  needed  to  determine  the  point  of  breaking 
and  the  wave  heights  after  breaking.  Though  a  number  of  formulas  for  doing 
this  have  been  developed,  there  is  not  as  yet  one  which  is  universally  applicabl 
or  accepted.  The  choice  of  a  breaking  criteria,  although  somewhat  arbitrary, 
must  be  made  with  care  since  it  determines  the  width  of  the  surf  zone  and  thus 
controls  the  set-up .  The  simplest  breaking  criteria  is  that  predicted  by 
solitary  wave  theory. 

(~)  =  constant  =  .78  (2.38) 

b 

where  the  subscript,  b,  denotes  the  value  at  breaking.  There  is,  however, 
considerable  evidence  (Weggel,  1972)  that  this  is  an  oversimpli f ica t ion . 

Noda  al .  (1974)  used  a  modified  version  of  the  Miche  formula 

(?)  =  .12  tanh  (?)  (2.39) 

L  b  b 

both  to  predict  the  point  of  breaking  and  the  decay  of  the  wave  after  breaking. 
This  was  done  by  calculating  both  a  wave  height  from  Eq,  (2.37)  and  a  breaking 
height  from  Eq.  (2,39).  When  the  point  was  reached  where  the  wave  height  was 
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equal  to  or  greater  than  the  breaking  height,  the  wave  was  considered  to 
have  broken  and  the  wave  height  from  Eq .  (2.39)  was  used.  Eq .  (2.39)  is 
used  in  both  models  to  determine  broken  wave  heights. 

2 . 8  Lateral  Mixing 


The  non-linear  model  of  Ebersole  and  Dalrvmple  (1979)  includes  the 
effect  of  the  lateral  shear  stress  terms 

_  St.  t  St  , 

-D  £  ~_D  t 

p  3y  5  p  3x 

in  the  x  and  v  momentum  equations  (2.12)  and  (2.13)  respectively.  The  need 
for  these  terms  is  pointed  out  by  Longuet-Higgins  (1970a)  in  his  treatment 
of  the  longshore  currents  due  to  obliquely  incident  waves  on  a  plane  beach. 
Neglecting  the  effects  of  lateral  shear  stresses  led  to  nrediction  of  a 
longshore  current  distribution  with  a  discontinuity  at  the  breaker  line  and 
no  current  outside  the  surf  zone.  However,  physical  observation  in  both 
the  laboratory  and  field  indicate  that  mean  longshore  flows  are  present 
beyond  the  breaker  line.  Longuet-Higgins  (1970b)  presented  a  formulation 
which  included  the  effect  of  lateral  mixing  as  the  means  for  handling  the 
effect  of  lateral  shear;  the  shear  stresses  are  thus  based  on  turbulent 
Reynolds  stresses  proportional  to  the  local  gradient  of  the  mean  velocity. 
The  resulting  velocity  distributions  have  no  discontinuity  at  the  breaker 
line,  and  the  peak  velocities  are  shifted  shoreward  of  the  breaker  line. 
Figure  2-2  shows  a  comparison  of  a  theoretical  velocity  profile  to 
one  without  lateral  mixing  included. 
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Figure  2-2  Longuet-Higgins ’  Analytical  Solution  for  Oblique 
Wave  Attack  on  a  Plane  Bead) 


Following  the  derivation  of  Longue t-Higgins  (1970b),  the  lateral 
shear  stress  is  written  as 


=  -  0  ( c 


9iJ 

3v 


av 

x  Dx 


) 


(2. AO) 


Longue t-Higgins  argued  that  the  mixing  coefficient  must  tend  to 
zero  as  the  shoreline  was  approached  since  the  turbulent  eddies  responsible 
for  mixing  cannot  be  of  a  scale  greater  than  the  distance  to  shore.  He  assumed 
that  is  proportional  to  the  distance  offshore,  x,  multiplied  by  some  velocity 
scale  which  he  chose  to  be  /gh,  the  speed  of  a  wave  in  shallow  water  where  h 
is  the  local  water  depth.  Therefore,  can  be  written  as 

c  =  Nx  /gh 
x 


where  N  is  a  dimensionless  constant  whose  limits  Longuet-Higgins  gave  as 


0  £  N  1  0.016 
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In  this  model  the  coefficient,  c  ,  was  allowed  to  vary  linearly  with  x  to 

x 

some  value  around  the  breaker  line.  From  this  point  seaward  the  coefficient 
remained  at  this  constant  value.  The  reason  for  this  approxima t ion  is  that 
physically  there  must  be  some  limit  on  the  scale  of  these  eddies.  This  limit 
is  at  present  not  known.  The  coefficient,  ,  was  chosen  to  be  constant. 

The  values  of  N  and  are  chosen  during  the  calibration  of  the  model. 

2 . 9  Wave  Height  Decay 

In  all  applications  of  the  circulation  model  to  date,  calculation 
of  wave  parameters  has  been  confined  to  a  region  within  several  wavelengths 
of  the  shore.  Over  this  distance,  wave  energy  decay  due  to  interaction  with 
the  bottom  is  not  likely  to  be  significant,  except  possibly  in  the  case  of 
waves  propagating  over  a  soft  mud  bottom  (Dalrvmple  and  Liu  (1978)),  which 
is  not  treated  here.  However,  the  circulation  model  could  reasonably  be 
used  to  model  propagation  over  much  longer  offshore  extents  of  shallow  coastal 
waters.  At  some  point,  the  accumulated  effects  of  bottom  interaction  would 
result  in  wave  height  reductions  of  a  significant  degree.  In  order  to 
accurately  model  the  amount  of  wave  energy  available  for  driving  currents 
and  maintaining  mean  water  level  variations,  it  is  necessary  to  include  the 
effects  of  various  dissipation  mechanisms  in  the  equation  for  calculating 
wave  height  (2.37). 

Ebersole  and  Dalrymple  (1979)  included  in  the  nonlinear  computer  model. 


but  did  not  describe,  the  option  to  calculate  wave  energy  dissipation  due  to 
viscous  shear  in  the  bottom  boundary  layer,  and  due  to  the  effect  of  "pumping" 
of  water  through  the  permeable  sand  bottom  due  to  wave  induced  pressure 


gradients  at  the  water-bottom  interface.  At  present,  the  option  of  calculating 
wave  energy  decay  due  to  both  mechanisms  is  included  in  both  the  linear  and 
nonlinear  model. 


The  application  of  wave  damping  to  the  nearshore  circulation  model 
lias  been  described  by  McDougal  (1979),  who  included  the  energy  dissipation 
rate  z  in  Eq .  (2.14)  in  a  sediment  transport  model  based  on  the  linear  wave- 
current  interaction  model.  The  derivation  of  e  is  based  on  Liu  (1973). 

Let  z  be  given  by 
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P 


f  . 
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(2.41) 


where  c  is  the  dissipation 
P 

l ^  is  the  dissipation  rate 
treated  separately. 
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Liu  (1973)  solved  the  linear  problem  for  waves  over  a  porous  bed  of 

H 

infinite  depth.  If  the  wave  amplitude  a(x,y,t)  =  — (x,y,t)  is  assumed  to  be 
a  slowly  varying  function  of  the  form 

a(x,y,t)  =  aQe 


Liu's  solution  leads  to  the  following  form  for  a; 
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where  Q  is  given  by 


(2.42) 
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(2.43) 


and 


-10  2 

=  bed  permeability  *  10  M 
Pq  =  bed  porosity  '  0.6 
v  =  kinematic  viscosity  * 

For  these  values  of  K  ,  P  ,  and  v,  the  viscositv  term  in  Q  is  dominant,  and 

p  o  '  * 

a  may  be  reduced  to  the  form 


2  cosh  kh 
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(2.44) 


The  dissipation  rate  a  is  related  to  c  by 
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"p  Jt 
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(2.45) 


by 


The  rate  of  energy  dissipation  due  to  bottom  friction,  e^,  is  given 


ef =  (Tb  ’  u)\ 


(2  .46  ) 


where  is  the  bed  surface  area. 


Substituting  for  and  U  using  results  derived  in  the  absence  of  a 
mean  current,  we  obtain  the  form 

-  Pf  t  \ 2 

t  0  77  m 


where  is  the  maximum  wave  induced  velocity  at  the  bottom,  given  by  (2.26), 
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After  some  manipulation,  we  obtain  the  form 


_k 

6-r 


a  f  H 


(cosh  kh-cosh  kli) 


(2.47) 


The  total  energy  dissipation  including  the  effects  of  the  porous 
bottom  and  friction  is  given  by 


e  +  £_ 

p  f 


+  ^  f  H _ 

cosh^kh  ^  6  ’t  (cosh^kh-cosh  kh)  j 


t  •  E 


(2.48) 


The  option  of  including  these  effects  is  available  in  both  forms  of 
the  numerical  model  presented  here. 


3.  AN  OVERVIEW  OF  THE  LINEAR  AND  NONLINEAR  MODELS 

While  both  of  the  numerical  models  described  in  Chapter  III  are  based 
on  the  theory  outlined  in  the  previous  section,  each  model  contains  a  somewhat 
different  subset  of  the  overall  development .  In  this  section,  we  review  the 
differences  and  similarities  between  the  models  before  going  on  to  their 
numerical  formulation. 

The  intent  of  both  of  the  models  is  to  solve  for  the  mean  values  U, 

V,  and  p  at  each  grid  point  by  solving  the  time  and  depth  averaged  equations 
of  continuity  (2.11)  and  momentum  (2.12-13).  Each  model  utilizes  the  refrac¬ 
tion  scheme  of  Noda  et_  aj^ .  (1974),  as  represented  by  (2.35-37),  to  solve  for 
wave  angle  and  wave  height.  The  model  of  Birkemeier  and  Dalrvmple  (1976) 
treats  linearized  forms  of  Eqs.  (2.12-13),  obtained  by  dropping  the  convective 
acceleration  terms 
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£  (ij2d)  *  37  (UVD) 

from  t ho  x-momentum  equation*  and  the  terms 

T-  (UVD)  ;  ~  (V2D) 

from  the  v -momentum  equation.  The  model  of  Ebersole  and  Dalrvmple  (1979) 
retains  the  full  nonlinear  form  of  the  momentum  equations,  leading  to  the 
principle  theoretical  difference  between  the  models  as  we 1 1  as  explaining 
the  distinction  expressed  by  their  names.  The  mathematical  differences  in 
the  governing  equations  also  lead  to  the  requirement  of  significantly 
different  numerical  schemes,  discussed  in  Chapter  III. 

The  momentum  equations  contain  various  forcing  terms  which  are 
formulated  as  stresses  or  stress  gradients;  these  include  radiation,  surface 
and  bottom  stresses  and  a  lateral  stress  representing  the  effect  of  turbulent 
mixing.  The  models  treat  radiation  stresses  and  surface  wind  stresses  iden¬ 
tically.  The  linear  model  treats  bottom  stress  according  to  a  '‘weak  current" 
formulation  developed  by  LeBlond  and  Tang  (1974),  based  on  the  assumption 
that  the  stress  develops  principally  in  response  to  the  wave  orbital  motion. 
The  nonlinear  model  uses  a  more  exact  bottom  stress  formulation,  making  no 
assumption  as  to  the  relative  magnitude  of  wave  orbital  and  mean  current 
velocities.  This  distinction  between  the  models  is  more  historical  than 
essential;  the  linear  model  can  be  updated  to  include  the  more  exact  relation 
given  by  Eqs.  (2.24)  and  (2.25).  This  modification  has  been  used  by  Allender 
et  al .  (1981)  in  a  version  of  the  linear  model  described  here.  A  further 
possibility  would  be  to  represent  the  bottom  friction  using  a  "large-current" 
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formulation  developed  hv  hiu  and  DalrympU-  (IV78);  tin's  extension  has  not 


boon  i  nvvs  t  i  gat vd  . 

Hie  remaining  d  i  st  i  :ic  t  ion  between  the  models  stems  from  the  neglect 
of  the  "lateral  mixing1*  terns  Kqs.  (2.12-13) 

-u  * \  m  -D  'lZ 

' V  '  e 

in  the  linear  model,  and  their  inclusion  in  the  nonlinear  model.  This 
difference  is  again  historical  rather  than  essential;  these  terms  could  be 
included  in  a  modified  version  of  the  linear  model,  although  their  inclusion 
would  require  more  effort  than  the  bottom  friction  modification  discussed 
above.  No  investigation  of  the  effect  of  including  lateral  mixing  in  the 
linear  model  has  been  made  to  date. 
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Chapter  III 


FORMULATION  OF  THE  NUMERICAL  MODELS 


1.  INTRODUCTION 

In  the  previous  chapter  a  mathemat ical  formulation  of  t\w  governing 
equations  and  associated  boundary  conditions  for  the  problem  of  nearshore 
wave- induced  circulation  has  been  outlined.  In  this  chapter,  tile  numerical 
formulations  of  the  mathematical  model  are  reviewed,  corresponding  to  the 
work  of  Birkemeier  and  Dalrymple  (1976)  and  Lbersole  and  Dalrvmple  (1979). 

The  two  models  reviewed  shall  be  referred  to  as  the  "linear  model" 
and  the  "nonlinear  model."  The  models  contain  significant  differences  as 
well  as  significant  similarities  in  their  numerical  formulation  as  well  as 
in  their  underlying  mathematical  formulation.  Differences  in  the  models 
will  be  discussed  below  following  a  review  of  the  overall  structure  common 
to  both  models. 

Both  models  reviewed  here  are  finite-difference  approximations  to 
a  set  of  three  first  order  hyperbolic  differential  equations  consisting  of 
the  continuity  equation  (2.11)  and  the  x  and  y  direction  momentum  equations 
(2.12,  2.13),  with  associated  unknowns  U(x,v,t),  V(x,v,t)  and  n(xtv,t).  The 
quantities  of  wave  height  H(x,v,t)  and  wave  angle  f(x,y,t)  are  solved  lor 
using  the  refraction  scheme  Eq .  (2.35)  and  the  wave  onergv  equation  (2.37). 
models  as  developed  contain  the  flexibility  to  handle  unsteady  response  to 


The 
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time  varying  incident  wave  conditions;  however,  vt-  will  be  concerned 
exclusively  with  representing  the  response  to  a  constant  wave  climate. 

In  this  guise,  the  models  represent  iterative  schemes  to  determine  the 
response  to  a  steady  state  exciting  force,  and  updated  unknowns  may  be 
regarded  as  iterated  values  rather  than  advanced* in -time  values. 

Listings  of  the  computer  programs  for  the  linear  and  nonlinear  models 
are  given  in  Appendix  I  together  with  some  notes  on  running  the  programs. 

2.  THE  GRID  SCHEME  AND  LOCATION  OF  THE  UNKNOWN  QUANTITIES 

The  first  step  in  constructing  a  finite  difference  model  lies  in 
choosing  a  network  of  discrete  grids  overlaying  the  physical  domain  of 
interest.  The  grid  used  by  the  present  models  is  that  of  Noda  et  al .  (1974), 
as  illustrated  in  Figure  3-1.  The  physical  domain  is  divided  into  regular 
grids  of  longshore  extent  Ay  and  offshore  extent  Ax.  The  topography  is  assumed 
to  be  periodic  in  the  longshore  v  direction  with  a  length  \  given  by 

A  =  (N  *  1)  •  Ay 

Various  requirements  affecting  the  choice  of  grid  are: 

1.  The  grid  must  extend  offshore  far  enough  to  remove  the  offshore 
region  of  the  domain  from  the  influence  of  currents  driven  by 
the  surf  zone,  and  to  allow  for  the  specification  of  a  uniform 
longshore  depth  which  will  not  significantly  alter  the  wave 
refraction  results  In  the  nearshore. 

2.  The  grid  mesh  must  be  fine  enough  to  resolve  the  surf  zone  in 
a  reasonable  manner. 
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3.  In  the  event  that  a  single  physical  feature  isolated  in  the 


longshore  direction  is  to  be  modelled,  the  longshore  extent 
of  the  grid  must  be  sufficiently  large  to  isolate  the  physical 
system  from  the  effect  of  images  created  by  the  longshore 
periodicity  requirement. 

In  the  event  that  large,  non-per iodic  physical  features  creating 
significant  current  patterns  seaward  of  the  breaker  line  are  present,  require¬ 
ments  (1)  and  (3)  lead  to  the  choice  of  a  large  grid,  whereas  require¬ 
ment  (2)  can  lead  to  choice  of  quite  small  offshore  grid  spacings  on  steep 
beaches.  The  resulting  grid  will  often  contain  a  large  number  of  elements, 
leading  to  the  requirement  of  large  amounts  of  computer  time  to  solve  the 
complex  set  of  governing  equations. 


The  fixed  and  variable  quantities  are  defined  in  relation  to  the 


grid  structure  as  shown  in  Figure  3-2.  The  quantity  ^  is  defined  at 
the  grid  center  and  represents  any  of  the  set  of  quantities 


A.  , 


=  (H, 


k,  n,  s,  D, 


V 


T  }  .  . 

S  1,J 


(3.1) 


Velocities  U  and  V  are  given  at  grid  sides.  This  formulation  relates 
naturally  to  the  conservation  law  scheme  of  relating  changes  of  a  quantity 
in  a  bounded  region  to  the  fluxes  of  that  quantity  across  the  bounding 
surface,  and  is  therefore  physically  correct  in  an  heuristic  sense,  as 
well  as  possessing  whatever  degree  of  numerical  accuracy  consistent  with 
the  chosen  difference  schemes. 
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OFFSHORE 


Figure  3-2 


Velocity  Components  for  Grid  Block  (i,j). 
Other  Variables  (D,  Sxx>  S  ,  H,  etc.)  Are 
Determined  at  Grid  Center. 


3.  BOUNDARY  CONDITIONS 


In  order  to  formulate  the  numerical  models,  a  scheme  was  established 
which  incorporated  the  lateral  periodicity  requirements  mentioned  above  and 
the  no-flow  requirements  at  the  shore  and  the  offshore  grid  row,  as  mentioned 
in  Chapter  II* 


The  user  defined  depth  grid  is  established  for  the  M  by  N 
portion  of  the  total  grid  shown  in  Figure  3.1,  with  the  requirement  that 


to  satisfy  the  periodicity  requirement. 

N  +  1  and  the  N  +  2  columns  according  to 


(3.2) 


The  grid  N  then  extended  to  the 


°i  ,N+1.  Di  ,  2 


D  ,  AT  ry  =  D  „ 

l , N+2  i , 3 


(3.3) 


Similar  periodicity  requirements  are  met  by  all  calculated  quantities  (the 

A.,  in  Eq.  (3.1)).  The  models  calculate  values  of  the  A.,  quantities  from  the 
1J  ij 

j  =  2  row  to  the  j  =  N  row.  Periodicity  is  then  established  by  setting 


A.  *  =  A,  „ 
1,1  i,N 

A  =  A 

i,N+l  i,2 

A  =  A 

i ,  N-f  2  i,3 


(3.4) 


The  velocities  U  and  V  are  treated  somewhat  differently  in  that  the  calculation 


are  performed  from  the  j  =  3  to  the  j  =  N+l  grid  rows,  with  periodicity 


established  accordingly. 

At  the  offshore  grid  row  (i  =  M) ,  a  no-flow  condition  is  applied 
for  U  and  V.  This  condition  serves  to  bound  the  calculated  velocities  during 
initial  start-up  of  the  model;  however,  the  condition  is  artificial  in  that 
the  offshore  side  of  the  modelled  area  becomes  a  solid  vertical  boundary. 

This  treatment  leads  to  potential  seiching  in  the  model.  The  presence  of  a 
seiching  mode  in  the  model  calculations  has  been  discussed  thoroughly  by 
Birkemeier  and  Dalrymple  (1971)  and  Ebersole  and  Dalrymple  (1979),  who 
have  shown  that  the  period  of  seiching  can  be  accurately  predicted  by  the 
shallow  water  formulas  for  known  topographies.  For  example,  for  a  basin  of 
triangular  cross-section,  Wilson  (1966)  predicts  the  period  of  the  first 
mode  oscillation  by 


where 


T 


3.28L 


(3.5) 


T  =  period  of  oscillation  in  the  basin 

L  =  length  of  the  basin  =  (M  -  1)  •  Ax 

h  ==  maximum  still  water  depth  in  the  basin, 

max 


Similarly,  a  no  flow  condition  for  the  shore-normal  velocity  U  is 
applied  at  the  onshore  grid  row,  whose  position  remains  fixed  during  a  model 
run.  The  linear  model  originally  described  by  Birkemeier  and  Dalrymple  (1976) 
included  a  provision  for  flooding  shoreward  dry  grids  in  order  to  model  the 
effect  of  wave  set-up;  model  results  have  been  seen  to  be  somewhat  insensitive 
to  the  correction  provided  by  this  provision.  Flooding  is  not  included  in  the 
nonlinear  model  of  Ebersole  and  Dalrymple  or  in  the  present  version  of  the 


linear  model. 
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In  addition  to  the  lateral  boundary  conditions,  the  system  of  hyperbolic 
equations  requires  initial  conditions  for  a  complete  solution.  The  models  are 
started  from  a  state  of  rest  with  no  wave  field  present.  In  order  to  reduce  the 
effect  of  seichirig,  the  wave  height  H  at  the  offshore  grid  is  brought  up  to  its 
full  value  gradually  according  to 

H  =  H  tanh 
o 

where 

t  =  model  time 

T  =  arbitrary  fixed  time  period. 

Wave  height  is  typically  allowed  to  build  up  for  400  model  iterations.  It  is 
also  noted  that,  in  principal,  the  seiching  effect  could  be  removed  by  setting 
n  to  zero  at  the  offshore  boundary,  rather  than  U. 


(3.6) 


4.  FINITE  DIFFERENCE  OPERATORS  AND  NOTATIONS 

The  derivations  of  the  numerical  models  given  by  Birkemeier  and 
Dalrvmple  (1976)  and  Ebersole  and  Dalrymple  (1979)  differ  significantly  in 
notation.  For  this  reason,  a  standardized  notation  is  presented  here.  We 
retain  as  far  as  possible  the  terminology  of  Ebersole  and  Dalrymple  (1979) 
in  describing  both  the  linear  and  nonlinear  models. 

The  numerical  formulations  require  both  averaging  and  differencing 
operations.  Let  F(x,y,t)  be  an  arbitrary  function  varying  in  space  and  time. 
The  average  of  the  function  over  one  grid  spacing  is  given  by 

F (x , y , t )  X  =  y  (F(x  +  -y  *  y,t)  +  F(x  -  ,y,t)l  (3.7) 

for  an  average  in  the  x-direct ion .  Successive  averaging  is  given  by 


F(x,y,t)  Xy  i£  F(x,y,t)  X 


(3.8) 


First  derivatives  can  be  given  in  terms  of  forward  and  backward  differences 
over  a  single  grid  space,  or  as  central  differences  over  two  grid  spaces.  We 
define  the  forward  and  backward  difference  operators,  respectively,  as 

6x(F(x,y,t)}  =  ~  { F (x  +  Ax,y, t)  -  F(x,y,t)}  (3.9) 

UF(x,y,t)}  =  ~  {F(x,y,t)  -  F(x  -  Ax.y.t)  }  (3.10) 

and  the  central  difference  operator  as 

•MF(x,y,  t) }  =  2^  (F(x  +  Ax,y,t)  -  F(x  -  Ax,y,t)}  (3.11) 

We  also  define  an  auxilary  operator  which  represents  the  central  difference 
at  a  grid  center  based  on  values  at  the  grid  sides 

<5~{F(x,y,t)}=  — ;  F (x  +  ,y,t)  -  F(x  -  ^  ,y,t)}  (3.12) 

It  is  easily  shown  by  direct  substitution  that  the  6~  operator  is  related  to 
6^  through  the  relation 

6~{F(x,y,t)  X}  -  6-{F(x,y,t)}  (3.13) 

x  x 

The  nonlinear  model  makes  extensive  use  of  the  <S^  operator ,  while  the  linear 
model  is  formulated  more  conventionally  in  terms  of  the  standard  forward  and 
backward  differences  given  by  Eqs.  (3.9)  and  (3.10). 


5.  FINITE  DIFFERENCE  FORMS  OF  THE  GOVERNING  EQUATIONS  -  LINEAR  MODEL 

Before  applying  the  finite  difference  scheme  to  the  linear  formulation, 
we  rewrite  the  equations  here  for  clarity.  The  equations  of  continuity  and 
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x-  and  y-momentum  are  given  respectively  by 
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(3.14) 


(3.15) 


(3.16) 


In  deriving  this  set  of  equations,  we  have  tacitly  assumed  that 
variations  of  n  with  respect  to  time  are  small  in  comparison  to  variations 
in  U  and  V.  As  we  are  striving  for  a  steady  state  solution,  the  model  results 
should  not  be  sensitive  to  this  assumption.  Care  should  be  taken,  however,  in 
removing  time  derivatives  of  D  in  cases  where  time  dependent  results  are 
desired . 


In  keeping  with  standard  techniques  for  treating  first  order  linear 
hyperbolic  equations,  difference  equations  are  formulated  using  a  forward  time- 
forward  space  (FTFS)  scheme  for  the  continuity  equation  (3.14),  and  forward 
time-backward  space  (FTBS)  formulation  for  the  momentum  equations  (3.15)  and 
(3.16).  Recalling  that  depths  and  set-up  are  known  at  grid  centers,  while 
velocities  are  given  at  grid  boundaries,  equations  (3.14)  -  (3.16)  are  written 
in  finite  difference  form  as: 
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Note  that  ^(F)  yields  the  forward  difference 


6t(F)  = 


Fk+1  -  Fk 
At 


(3.20) 


where  k  and  k+1  are  successive  time  levels.  Expanded  forms  of  the  equations 
(3.17)  -  (3.19)  can  be  found  in  Birkemeier  and  Dalrymple  (1976).  At  each 
iteration,  the  values  of  U  and  V  are  updated  to  time  level  k+1  using  n  and 
forcing  terms  evaluated  at  time  level  k.  Then  n  is  updated  using  the  newly 
evaluated  values  of  U  and  V  at  time  level  k+1.  The  method  thus  requires  only 
one  time  level  of  storage  lor  each  of  the  unknown  quantities. 


Difference  forms  of  the  type  used  here  have  been  discussed  extensively 
by  a  number  of  authors  (see,  for  example,  Roache  (1976)).  The  explicit  method 
is  subject  to  certain  stability  conditions.  The  problem  of  obtaining  stability 
criteria  for  equations  with  non-constant  coefficients  is  in  general  unsolved 
to  date;  however,  we  can  make  a  reasonable  guess  based  on  constant-coefficient 
forms  of  the  governing  equations.  If  we  drop  the  forcing  terms,  assume 
D  '  hf  and  hold  h  constant,  we  obtain 


h  M  +  iX 

:>t 

"  h  Ux  +  3y; 

w  _ 

*t 

8  3x 

JV  _ 

*n 

it 

-8^ 

(3.21) 
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Cross-differentiating  to  eliminate  U  and  V,  we  obtain  the  second-order 
hyperbolic  equation  for  n; 


> 


•t 


(3.22) 


The  stability  criterion  corresponding  to  this  reduced  equation  is  given  by  the 
Courant  condition: 


't  _  ,  (3.23) 

‘2gh 

assuming  that  v  =  Ax.  Generalizing  to  the  full  model,  the  stability  criterion 
(3.23)  can  be  extended  to  include  the  effect  of  rectangular  grids  and  the 
advection  of  disturbances  by  the  mean  currents. 


9  9 

y(\x)^  +  (Ay)" 

v2^h  +  Jij| 


(3.24) 


The  stability  criterion  basically  states  that  time  steps  in  the  model  may  not 
be  so  large  as  to  allow  long  wave  disturbances  to  pass  completely  through  a 
grid  in  one  iteration.  Violation  of  the  criterion  leads  to  rapidly  growing 
instability  of  the  numerical  solution.  In  practice,  time  steps  on  the  order 
of  0.25  times  the  maximum  At  are  used. 


b.  FINITE  DIFFERENCE  FORMS  OF  THE  GOVERNING  EQUATIONS  -  NONLINEAR  MODEL 

The  set  of  nonlinear  momentum  equations  (2.12)  -  (2,13),  together  with 
the  continuity  equation  (2.11),  require  a  more  careful  development  in  difference 
form  than  the  previously  described  linear  model.  Nonlinear  formulations  are 
in  a  practical  sense  subject  to  a  number  of  instability  mechanisms  which  are 
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not  strictly  related  to  stability  criteria  for  the  corresponding  linear 
formulations,  as  discussed  by  Roache  (1976). 

The  nonlinear  model  has  been  formulated  using  a  method  which  has  been 
applied  successfully  to  geophysical  and  tidal  models  by  Lilly  (1965)  and 
Blumberg  (1977).  Using  the  difference  and  average  operators  defined  in 
section  4,  the  equat ions (2 . 11-13)  are  rewritten  as 

Cent i nu i tv : 

6;{"n  C1  +  <5~{D  Xl0  +  6~{D  yV}  =  0  (3.25) 

t  x  y 

X-  Mo  men  turn 

— -  t  ’3 — 7-  x  _  —  -  --  x  _ 

6 ' { D  XU  }  +  ^'{D  XU  U  X}  +  6- { D  >V  U  y)  = 

t  X  V 


=  -gDX6{n)+—  (t  x  -  t  -  6  ~ { s  X' }  - 

°  x  i'  sx  bx  v  Jxv 


-  6-{S  } }  +  D  X  &~it  5~1U)  +  e  Xy  6-{V} 

x  xx  y  y  y  x  x 


(3.26) 


Y-Momen turn 


z —  t  z —  y__  Z — T~  -v— 

"(D  VV  }  +  6~{D  XU  V  y }  +  6'{D  yV  V  y) 

t  X  V 


__  —  i  f - y  —  y 

g  D  ^  6  in)  +  ~  U  -  x  -  6~{S  } 

*  y  P  V  sy  by  v  yy 


-  ~  {  S  xy}i  +  D  y  6~{ e  6  '  (U  }  +  c  xy  5 '  { V }  } 

x  xv  |  x  y  y  x  x 


(3.27) 


The  formula  for  lateral  mixing  given  by  Eq .  (2.46)  lias  been  substituted 
into  Eqs .  (3.26-27).  Note  that,  following  Eq .  (3.13),  the  central  difference  on 
the  time  averaged  wate  set-up  in  Eq .  (3.25)  leads  to 
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S~  In 


}  -  l n) 


—k+1  -k-1 

n.  .  -  n.  . 

-JjJ _ AjJ. 

2  At 


(3.28) 


The  governing  equations  in  this  case  then  contain  function  values  at  three 
times  levels,  given  bv  k+1,  k,  and  k-1.  The  equations  (3.25-27)  can  be  i ound 


written 

out 

with  respect 

to  the 

i > j  K<"id 

in  ] 

Ibersole  and  Da  1  rymp  1  e  (1979). 

These  three 

equations  can 

a  1  so 

be  written 

in 

the  following 

abbreviated  form. 

—k+l 
n .  . 

-k-1 

=  n.  .  + 

2  At 

h 

(3.29) 

Uk+1 

i » j 

=  A  uf-1 
h  J 

+  2. 

■  t  f2 

(3.30) 

vk+! 

=  B  Vk_1 

1 

+  2. 

•t  f3 

(3.31) 

where  A 

and 

B  are  functions  of 

the  depth 

alone  and  F^ ,  F..,, 

and  F^  are  functions 

of  all  the  variables  in  the  problem.  The  quantities  ,  F<;),  and  F^  contain 
quantities  evaluated  at  time  level  k  and  k-1. 


The  method  of  solution  for  Eqs.  (3.29-31)  is  based  on  a  leapfrog  scheme, 
and  thus  requires  the  storage  of  three  time  levels  of  values  for  all  the 
calculated  variables.  In  order  to  select  the  model,  a  single  step  is  taken 
based  on  a  forward  difference  formulation  similar  to  that  described  in  con- 


junction  with  the  linear  model. 

indicated  schematically  by 

— k+1 

The  forward 

— k 

difference 

formulation  can  be 

(3.32) 

ni»  j 

=  n.  .  + 
i » J 

At 

Fi 

uk+1 

=  cuk, 

hj 

+  A 

F4 

(3.33) 

vk+1 
i,  j 

=  D  V.  . 
hj 

+  A 

,t  f5  . 

(3.34) 

The  computat ional  method  is  schematized  in  Figure  3-3. 
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Figure  3-3.  General  Leapfrog  Solution  Scheme 

The  stability  criterion  for  the  present  scheme  was  assumed  to  be 
given  by  Eq .  (3.24),  which  has  physical  if  not  mathematical  significance  in 
the  present  situation.  In  practice,  as  is  the  case  of  the  linear  model,  it 
was  necessary  to  reduce  the  maximum  allowable  time  step  to  a  value 
significantly  below  the  one  predicted  by  the  stability  criterion. 

The  nonlinear  model  was  also  subject  to  a  second  computational  stability 
problem.  In  general,  computational  schemes  for  first  order  equations  which  are 
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centered  in  space*  and  time*  (CSCT)  exhibit  forms  of  unstable  behavior.  In 
particular,  in  the  context  of  the  first  order  parabolic  diffusion  equation, 
such  explicit  CISC F  schemes  can  be  shown  to  be  unconditionally  unstable.  In 
tile  present  case,  as  the  model  approached  a  steady  state,  the  solution  di¬ 
verged  into  two  disjoint  solutions;  one  associated  with  the  even  time  steps 
and  the  other  the  odd  steps.  These  solutions  oscillated  with  growing  ampli¬ 
tudes  about  the  steady  state  solution.  In  Roache  (1976),  the  author  referred 
to  this  as  time  splitting. 

To  alLeviate  the  problem,  a  leapfrog-backward  correction  scheme, 

Kurihara  (1965),  was  initiated  every  tenth  time  step.  The  scheme  is  shown 
below  as. 

It*  =  hk_1  +  2.t  Gk  (3.35) 

k  +  1  k 

h  =  h  +  .It  G*  (3.36) 

where  h  may  be  U,  V  or  n.  Equation  (3.35)  is  simply  the  leapfrog  calculation 
done  by  the  Eqs.  (3.29-31)  where  *  denotes  the  new  or  "interim''  time  level. 

Using  the  new  values  U,  V,  n  at  time  *,  the  functions  G*,  like  the  functions 
F^,  F?,  and  F^  from  Eqs.  (3.29-31)  are  calculated  and  used  in  Eq . 

(3.36)  which  is  merely  a  backwards  difference  in  time  to  the  level  k. 

This  scheme  was  chosen  because  it  damps  the  computational  mode  of 
the  solution  while  leaving  the  physical  mode  relatively  unaffected.  For  a 
more  in-depth  discussion  the  reader  is  referred  to  the  work  by  Kurihara.  With 
this  correction  scheme  included,  which  essentially  ’’ties”  the  solutions  together 
every  tenth  iteration,  the  model  proceeded  to  reach  a  steady  state  without  any 
further  time-splitting  instability. 
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7.  L'HE  NUMERICAL  SCHEME  FOR  RK FRACTION  AND  THE  WAVE  HEIGHT  FIELD 

The  *q nations  governing  wave  refraction,  wave  height,  and  wave  number 

are  given  by  Eqs.  (2.35),  (2,37),  and  (2.36)  respectively.  The  method  of 

solving  these  equations  for  the  unknown  d,  . ,  H .  .  and  k.  .  is  identical  in 

iij 

both  the  linear  and  nonlinear  models.  The  method  was  given  originally  by 
Noda  et  al.  (1974). 


If  Eq .  (2.36)  is  differentiated  with  respect  to  x  to  get 

3k 


3  k 
3x 


and  with  respect  to  y  to  get  — ,  these  quantities  can  be  substituted  into 
Eq .  (2.35),  which  can  then  be  written  in  the  following  form: 


A  3x  B  3v 


=  C 


(3.37) 


where  A,  B  and  C  are  functions  of  the  quantities  sin  6,  cos  6,  k,  h,  u  and  v. 

By  taking  a  forward  difference  in  x  to  approximate  —  and  a  backwards  difference 
in  v  to  approximate  Eq.  (3.37)  becomes: 

a  V 


i-j 


D+EQ.  .  .  -  F  0.^  . 

i,J-l  l+l, J 


(3.38) 


where  D,  E  and  F  are  now  iunctions  of  the  quantities  sin  0.  cos  9. 

i.J  hJ 


k.  . ,  h.  .  ,  u,  . ,  v.  ..  To  evaluate  sin  0.  .  and  cos 
i*J  i»J  i.J  i»J  itJ 


.  Noda  used  a  first 


order  Taylor  series  expansion  to  the  four  neighboring  grids  (i+l,j),  (i-l,j), 
(i,j+l)  and  (i,j-l),  summed  the  results  and  took  an  average  value. 


The  theta  field  0^  j  is  solved  for  in  the  following  way.  Snell’s 
Law  is  used  to  approximate  the  angles  at  the  offshore  row.  Working  shoreward 
Eq.  (3.38)  is  solved  for  in  a  row-by-row  relaxation  until  the  angles  converge 
to  their  correct  values  with  wave-current  interaction  included.  After  each 
updated  value  of  theta,  a  new  wave  number  must  be  solved  for. 
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Eq .  (2.36)  can  be  written  as 


E(k)  ::  {gk  tanh(kh)}1^  +  ukcosQ  +  vksinO  -  =  0  .  (3.39) 

To  solve  for  the  wave  number,  k,  after  each  updated  angle  is  found,  the  Newton- 
Raphson  Method,  or  "method  of  tangents",  is  used.  This  method  states  that 

E(k  ) 

k  =  k  ,  .  -  -f  — T 
new  old  E  (k  ,  , ) 
old 

Dif ferentiat ing  Eq.  (3.39),  k  is  iteratively  solved  for  until 

new  J 

|  k  -  k  .  .1  i  . 001 1 k  I  . 
new  old'  new 


The  wave  height  field  is  calculated  in  much  the  same  way  as  the  wave 

H  3R 

angle  field.  Multiplying  Eq .  (2.37)  by  --  and  letting  —  *  0,  the  energy 

Z  at 

equation  can  now  be  written  in  the  form 


A  —  +  B  =  C  H  (3.40) 

3x  3y 

where  A,  B  and  C  are  functions  of  u,  v  cos  0,  sin  G,  C  ,  Ax,  Ay  and  the  radiation 

S 

9H 

stresses.  Taking  a  forward  difference  in  x  to  approximate  —  and  a  backward 

3H 

difference  in  v  to  approximate  —  and  solving  for  H.  it  can  be  show  that 


H.  . 
hj 


D  H.  .  , 
i»J-l 


E  H .  . 


where  D  and  E  are  functions  of  the  same  quantities  as  A,  B  and  C.  Again  the 
offshore  row  of  wave  heights  are  obtained  from  shoaling  and  refraction  due  to 
Snell’s  Law  and  the  wave  height  field  is  determined  by  relaxing  row  by  row  in 
the  shoreward  direction.  On  each  row  a  solution  for  the  wave  height  is  reached 
when  |H  -  H  .,1  5  .01  H  .  After  each  updated  value  of  H.  a  breaking 
wave  height  is  also  calculated  from  the  breaking  criteria  given  by  the  Miche 
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formula 


If  the 
height 


(~r)  =  .12  tanh(kh), 

L  b 


calculated  H.  .  is  larger  than  the  allowable  breaking  height, 
1  *  1 

IL  j  was  set  equal  to  the  breaking  height. 


the 
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CHAPTER  IV 


CALIBRATION  OF  THE  NEARSHORE  CIRCULATION  MODELS 

In  order  to  make  the  circulation  model  more  generally  applicable 
to  prediction  of  field  conditions,  both  versions  of  the  model  were  cali¬ 
brated  using  field  data.  The  calibrations  consisted  of  determining  a 
monochromatic  deepwater  wave  condition  and  a  uniform  wind  condition,  and 
then  running  the  models  using  the  field  conditions  and  measured  bathymetry 
Adjustable  coefficients  were  then  tuned  to  obtain  the  best  possible  quali¬ 
tative  and  quantitative  fit  between  the  currents  predicted  by  the  model 
and  those  measured  in  the  field. 

In  this  chapter,  the  field  data  set  used  for  calibration  is 
described.  Results  for  the  circulation  models  are  then  shown  for  a  range 
of  coefficient  values. 

1.  FIELD  DATA  USED  FOR  CALIBRATION 

Field  data  used  during  calibration  of  the  nearshore  circulation 
models  was  obtained  from  the  results  of  the  Nearshore  Sediment  Transport 
Study  (NSTS)  experiment  conducted  at  Torrey  Pines  beach,  near  La  Jolla,  Ca 
(see  Figure  4-1)  in  November  -  December  1978.  The  results  of  this  study 
were  chosen  as  being  applicable  to  calibration  of  the  models  for 
several  reasons*  Torrey  Pines  beach  is  a  long,  straight  beach  with  fairly 
uniform  and  parallel  offshore  contours.  The  field  bathymetry  was  thus 
easily  adapted  into  the  model's  scheme  of  longshore  periodicity. 


49 


Figure  4 -1 . 


Location  of  NSTS  Torrey  Pines  Experiment . 
(i-iMiii  (miIiL  ,  1 1 1 7 9 ) 


The  experiment  itself  produced  a  detailed  picture  of  nearshore 
currents  and  waves,  with  up  to  22  current  meters,  10  pressure  gages,  and 
7  wave  staffs  being  operated  simultaneously.  Thus  the  resolution  needed 
to  accurately  fit  the  model  predictions  was  present  in  the  field  data. 
Gable  (1979)  gives  a  detailed  description  of  the  site,  instrumentation, 
and  conduct  of  the  NSTS  experiment.  The  arrangement  of  fixed  instruments 
used  in  the  experiment  is  shown  in  Figure  4 -2 .  Angles  and  distances  used 
in  this  report  will  be  with  respect  to  the  baseline  (0  offshore  distance) 
in  Figure  4-2.  Two  separate  bathymetry  maps  are  shown  in  Figures  4-3 
and  4-4.  It  is  noted  that  the  survey  of  Nov.  9,  1978  indicates  the  pres¬ 
ence  of  a  shallow  channel  oriented  almost  perpendicular  to  the  baseline 
at  about  40  m  left  of  the  main  range  (0  m  alongshore) .  This  feature  is 
not  present  in  the  Nov.  18,  1978  survey,  which,  on  the  whole,  exhibits 
greater  random  fluctuation  in  the  location  of  the  contours.  Both  the 
channel  in  the  Nov.  9  bathymetry  data  and  the  unevenness  in  the  Nov.  18 
data  appear  to  be  transient  features,  as  will  be  discussed  below. 

1.1  Choice  of  Field  Data  for  Calibration 


Several  requirements  were  chosen  in  order  to  determine  a  valid 
set  of  field  data  for  comparison  to  the  numerical  model.  First,  the 
numerical  models  in  their  basic  state  are  designed  to  be  run  with  a  mono¬ 
chromatic  deepwater  wave  condition.  It  was  therefore  required  that  the 
wave  field  for  the  chosen  data  be  narrow  banded  both  in  a  frequency  and 
directional  sense.  The  presence  of  a  second  wave  component  at  a  different 
direction  than  the  primary  component  introduces  a  forcing  mechanism 
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Figure  4-3.  Nearshore  Bathymetry,  Nov.  9,  1978 
(from  Gable  1979) 
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Figure  4-4.  Nearshore  Bathymetry,  Nov.  18,  1978 
(from  Cable  1979) 
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which  would  tend  to  produce  rip  cells  on  a  plane  beach,  as  shown  by 
Dalrymple  (1975).  Just  such  a  case  was  treated  by  a  simplified  refraction 
model  in  Ebersole  and  Dalrymple  (1979) ,  but  the  capability  to  handle  the 
condition  has  not  been  retained  in  the  final  model  versions  since  the 
refraction  scheme  used  does  not  model  the  interaction  of  different  waves. 

In  addition,  the  use  of  a  17  minute  average  of  the  field  measured  currents 
in  order  to  obtain  fairly  steady  conditions  precluded  the  modelling  of 
essentially  transient  circulation  phenomenon. 

Secondly,  since  it  was  desired  to  calculate  a  steady  state  wave 
and  current  field,  it  was  required  that  wave  and  wind  conditions  be  nearly 
steady  over  the  interval  of  averaged  field  measurements.  This  required 
either  a  quiet  or  unidirectional  wind  field  as  well  as  fairly  steady  wave 
direction,  height  and  period. 

The  field  data  set  available  from  the  NSTS  results  was  roughly 
screened  on  the  basis  of  the  beach  observations  given  in  Gable  (1979).  In 
particular,  observation  of  long  crested  waves  indicated  the  presence  of 
narrow  banded  directional  spectra.  It  was  found  that  the  NSTS  data  presented 
a  problem,  in  that  almost  all  data  sets  exhibited  at  least  some  degree  of 
bi-directionality,  with  waves  approaching  from  both  the  north  and  south. 

The  data  chosen  for  initial  calibration  has  a  somewhat  smaller  wave  component 
from  the  north  than  other  records,  but  is  still  suspect  as  a  valid  calibration 
standard. 

Data  Set  Number  1 

The  NSTS  data  tapes  divide  the  current  and  pressure  records  into 
17  minute  segments.  The  first  data  set  chosen  for  use  in  cailbration. 


and  subsequently  used  for  the  majority  of  calibration  runs,  was  taken  from 


the  second  17  minute  segment  of  the  test  of  Nov.  10,  1978.  The  17  minute 
segments  were  further  subdivided  into  four  4.25  minute  segments,  and  averaged 

current  fields  were  plotted  for  each  4,25  minute  segment.  Two  representative 
plots  for  the  first  data  set  are  shown  in  Figure  4 -5a  and  b. 

Based  on  an  average  of  the  4.25  minute  average  velocity  profiles, 
an  average  17  minute  velocity  profile  was  constructed  for  data  values  on  the 
main  range  (offshore  array  of  meters)  and  is  shown  in  Figure  4-6.  The  re¬ 
sulting  velocity  profile  indicates  a  longshore  current  of  about  8  cm/sec 
magnitude  offshore  of  the  surf  zone.  This  current  is  too  far  outside  the 
breaker  line  to  be  surf  driven,  and  may  be  due  to  the  presence  of  a  tidal 
or  seasonal  current.  The  effect  of  this  current  on  the  model  results  is 
discussed  below. 


set  No. 


The  monochromatic 

deepwater  wave  conditions  determined  for  data 

1  were  as  follows 

Wave  period 

T  : 

14.52  seconds 

Wave  height 

H 

o 

0.62  meters 

Wave  angle 

A  : 

165.7°  measured  clockwise  from  +x 

(offshore  direction). 

The  uniform  wind  conditions  were: 

Wind  speed  W  :  7.2  meters/second 

Wind  angle  a  :  294.5°  measured  clockwise  from  -x 

(onshore)  direction. 
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Data  Set  Number  2 


Data  Set  Number  1  represented  the  only  strongly  unidirectional, 
narrow  banded  wave  spectrum  contained  in  the  Torrev  Pines  data.  The  re¬ 
maining  data  sets,  even  in  the  absence  of  locally  generated  short  wind  waves, 
exhibited  at  least  a  strong  tendency  towards  bidirectionality,  with  spectral 
peaks  nearly  symmetrically  placed  about  the  shore-normal  axis.  Since  it  is 
anticipated  that  measured  wave  fields  in  general  would  seldom  exhibit  the 
narrow  banded,  unidirectionality  required  as  input  into  the  model,  a  data 
set  was  constructed  based  on  averaging  over  the  parameters  of  the  measured 
wave  field  in  order  to  test  the  response  of  the  model  using  a  "best  guess'1 
for  the  desired  input  parameters. 

The  data  chosen  for  this  test  was  from  the  eighth  17-minute  time 
period  of  the  November  15  data.  Wave  energy  was  contained  in  two  directional 
peaks  in  a  narrow  frequency  band  (Figure  4-7).  Using  directional  spectra 
supplied  by  Pawka  (1980),  the  data  set  was  constructed  by  summing  the  mean 
energy  density  over  all  frequencies  in  the  dominant  spectral  peak.  The 
wave  angle  was  chosen  as  the  average  of  the  dominant  angle  at  each  frequency 
weighted  by  the  energy  densities.  The  peak  energy  containing  frequency  was 
taken  as  the  dominant  frequency.  The  resulting  data  set  follows: 


Wave 

Period 

T 

:  14.71  seconds 

Wave 

Height 

H 

o 

:  0.43  meters 

Wave 

Angle 

A 

:  168°  measured  clockwise  from  -fx 

(offshore  direction) 
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Figure  4-7.  Directional  Distribution  of  Wave  Energy  at  Peak  Frequencies 
Nov.  1  5,  19  78  (supplied  by  S.  Pawka  ,  1980) 
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There  was  no  significant  steady  wind  component. 

Wind  speed  W  :  0.0  meters/ second 

Wind  angle  at  :  not  applicable. 

The  resulting  data  set  has  an  indeterminant  connection  with  the 
measured  field  data,  and  at  best  the  model  could  be  expected  to  exhibit 
results  in  qualitative  agreement  with  the  field  data. 

Figure  4-8  shows  a  plot  of  average  longshore  current  over  the  17 
minute  period  of  data  collection.  Individual  plots  (not  shown)  of  4.25 
minute  averages  show  a  great  deal  of  scatter  in  peak  velocities  and  shape 
of  the  profile,  indicating  some  unsteadiness  in  the  current  field,  as  would 
be  expected  due  to  the  bi-directional  wave  field. 

2.  CALIBRATION  OF  THE  MODELS 

Both  the  linear  and  nonlinear  models  have  a  bottom  friction  factor 
f  and  the  Van  Dorn  coefficients  and  as  adjustable  parameters.  In 
addition,  the  nonlinear  model  has  adjustable  coefficients  of  lateral  eddy 
mixing  in  both  the  longshore  and  offshore  directions. 

The  Van  Dorn  wind  stress  formulation  used  in  the  models  is  intended 
to  represent  the  transfer  of  momentum  from  a  wind  field  to  the  water  column, 
leading  to  a  wind-induced  longshore  current  for  wind  stress  applied  in  the 
long  shore  direction,  and  wind  set  up  for  wind  blowing  towards  shore.  Wind 
set  up  cannot  currently  be  accurately  calibrated  in  the  present  models  due 
to  the  artificial  barrier  at  the  offshore  grid  row;  no  additional  water  can 
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Offshore  Distance  from  NSTS  Baseline  (m) 


enter  the  model  after  the  start  of  the  ran  with  a  given  depth  grid.  The 
parameters  used  in  the  Van  Dorn  formula  have  been  tested  in  previous  large 
scale  models  (Pearce*  1972)*  and  have  been  found  to  be  satisfactory. 

With  the  elimination  of  the  wind  stress  coefficients,  only  the  bottom 
friction  coefficient,  f,  remained  to  be  calibrated  in  the  linear  model.  The 
procedure  used  was  to  choose  a  value  of  f  for  both  models  based  on  a  comparison 
of  field  data  with  the  linear  model.  The  value  of  f  chosen  was  then  used  in 
the  nonlinear  model  as  a  first  approximation,  and  values  of  the  eddy  mixing 
coefficients  were  adjusted  to  again  obtain  a  best  fit  between  the  results  of 
the  nonlinear  model  and  the  field  data.  It  should  be  noted,  however,  that 
there  is  no  _a  priori  reason  that  both  models  should  behave  optimally  with 
the  same  value  of  f. 

2. 1  Linear  Model  Calibration 

Runs  of  the  linear  model  were  conducted  using  the  deepwater  wave 
conditions  of  Data  Set  1  and  the  measured  bathymetry  of  Nov.  9,  1978.  For 
all  values  of  f  chosen,  rip  currents  formed  near  the  shallow  channel  seen  in 
the  Nov.  9  bathymetry.  These  rip  currents  were  not  apparent  in  the  Nov.  10 
currents  (Figure  4-5),  where  they  would  affect  the  main  range  velocity  profile. 
An  example  of  the  current  field  calculated  using  the  Nov.  9  bathymetry  is 
shown  in  Figure  (4-9).  The  value  of  f  for  the  run  shown  was  0.01. 

It  was  tentatively  concluded  at  this  point  that  the  shallow  channel 
seen  in  the  Nov.  9  data  was  a  transient  feature  and  was  not  present  on  Nov.  10, 
juding  from  the  uniformity  of  the  measured  currents.  It  was  decided  to  run  a 
model  using  a  bathymetry  based  on  a  longshore  average  of  the  Nov.  9  bathymetry 
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data.  The  longshore  extent  of  the  beach  was  arbitrarily  chosen  as  200 
meters.  The  profile  is  shown  in  Figure  4-10. 


Using  the  longshore  averaged  bathymetry,  model  runs  were  conducted 
for  a  range  of  f  values.  Velocity  profiles  for  values  of  f  of  0.01,  0.015 
and  0.02  are  shown  in  Figure  4-11,  in  comparison  to  the  current  distribution 
measured  in  the  field.  In  addition,  a  "corrected”  field  current  distribution 
constructed  by  subtracting  the  0.08  m/sec  offshore  current  is  shown  as  the 
dashed  line. 

Figure  4-11  shows  that,  by  using  a  value  of  f  equal  to  0.015,  the 
linear  model  closely  predicts  the  maximum  velocities  and  the  general  velocity 
distribution  in  the  surf  zone.  The  longshore  current  predicted  by  the  model 
dies  off  more  quickly  offshore  due  to  the  absence  of  lateral  mixing  effects 
in  the  linear  model. 

2 . 2  Nonlinear  Model  Calibration 

The  nonlinear  model  was  run  using  Data  Set  1  and  the  value  of 

f  =  0.015  obtained  from  the  linear  model  calibration.  Three  sets  of  mixing 

coefficients  e  (or  N)  and  £  were  used  (see  Equation  (2.40)). 
x  y 

N  £ 

0.000  0.00  no  mixing 


0.0025 


0.25 
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Figure  4- 1 1 


Variation  of  Longshore  Current  with  Friction 
Factor  in  the  Linear  Model. 


Results  for  f  =  0.015  are  shown  in  Figure  4-12.  It  was  found  that 
the  inclusion  of  mixing  satisfactorily  extended  the  velocity  profile  in  the 
offshore  direction,  but  in  both  cases,  including  mixing,  the  current  magni¬ 
tudes  in  the  surf  zone  were  underpredicted.  The  model  was  therefore  retested 
using  a  smaller  value  for  the  friction  coefficient,  f  =  0.01.  Results  are 
shown  in  Figure  4-13.  The  smaller  value  of  f  is  seen  to  correct  for  the 
underprediction  of  longshore  current.  Coefficient  values  for  the  nonlinear 
model  were  chosed  based  on  the  second  set  of  results, 

f  =  0.01 
N  =  0.0025 

e  *  0.25 

y 

2 , 3  Response  of  Both  Models  Using  Data  Set  2 

The  linear  and  nonlinear  models  were  run  using  Data  Set  2  and  the 
calibrated  coefficients  obtained  above.  Results  for  the  linear  model  are 
shown  in  Figure  4-14.  Nonlinear  model  results  were  similar  to  the  linear 
model  results.  Both  models  were  seen  to  overpredict  the  maximum  longshore 
current  in  comparison  to  the  averaged  field  data,  and  to  underpredict  the 
offshore  extent  of  the  longshore  current.  It  is  likely  that  a  better  fit 
to  the  field  data  could  be  obtained  by  artificially  increasing  the  lateral 
mixing  in  the  nonlinear  model.  However,  the  field  data  represents  the 
average  of  a  complex  flow  pattern,  where  the  effect  of  mixing  in  the  averaged 
data  is  induced  by  large  scale  phenomenon  not  likely  to  be  found  in  the  steady 
current  fields  induced  by  a  unidirectional  wave  field.  Since  it  is  likely 
that  the  magnitude  of  artificial  mixing  would  vary  greatly  as  a  function  of 
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Figure  4-12  Variation  of  Longshore  Current  with  Lateral 

Mixing  in  the  Nonlinear  Model, 
f  =  0.015 
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Variation  of  Longshore  Current  with  Tiateral 
Mixing  in  the  Nonlinear  Model, 
f  =  0.010 
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Figure  4-14 
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wave  directionality  and  spatial  complexity,  it  would  be  unjustified  to  alter 


the  parameters  of  the  model  to  fit  a  single  case.  It  may  become  possible  at 
some  future  date  to  estimate  mixing  parameters  based  on  the  charac teri st ics 
of  the  random  wave  field. 

A  plot  of  the  velocity  vectors  for  data  set  2  obtained  using  the 
calibrated  model  and  the  November  18  bathymetry  is  shown  in  Figure  4-15. 

As  in  the  November  10  case,  the  circulation  pattern  is  seen  to  be  sensitive 
to  the  bottom  variations.  This  result  is  probably  due  to  the  low  values 
obtained  for  the  friction  coefficient  f. 

3,  DISCUSSION  OF  THE  CALIBRATED  COEFFICIENTS 

The  value  of  the  friction  factor  f  obtained  in  this  study  differs 
significantly  from  the  value  used  in  previous  work  and  retained  by  Allender 
e^t  al_.  (1981).  Based  on  a  bottom  shear  stress  relation  given  by  Eqs.  (2.27- 
2.28),  Birkemeier  and  Dalrymple  (1976)  chose  a  value  of  C^  =  0.01,  which 
corresponds  to  a  choice  of  f  =  0.08.  This  value  of  f  is  carried  over  into 
the  results  discussed  in  the  next  chapter.  However,  model  calibration  has 
indicated  that  the  value  of  f  is  more  of  the  order  0.01-0.02,  with  values 
of  0.015  and  0.01  chosen  for  the  linear  and  nonlinear  models  respectively. 

It  is  felt  that  this  alteration  in  choice  of  the  friction  factor  requires 
some  discussion. 

Many  formulae  exist  to  calculate  the  bottom  friction  factor  under 
waves,  the  most  successful  being  those  of  Jonsson  (1966)  and  Kajiura  (1968). 
Writing  Kajiura* s  formula  in  a  form  given  by  Dalrymple  and  Lozano  (1978), 
we  obtain 
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cf  = 


=  1.41 


4nd 


T/  i  2.1/2 


2/3 


where  d  is  the  median  sand  grain  dimeter  and  k  is  the  breaking  index.  The 
value  of  f  at  the  breaker  line  represents  a  reasonable  choice  for  an  average 
value  over  the  surf  zone  region.  For  Torrey  Pines  Beach,  d  is  approximately 
given  by 


d  -  d_n  =  0.27  mm 

oO 

For  the  field  data  used,  is  taken  as  approximately  1.5  m .  This  yields  a 
value  of 


Cf  =  .0025  or  f  =  0.02 

The  values  of  f  obtained  during  calibration  are  thus  of  the  correct  order  of 
magnitude  for  flows  over  a  planar  bed  with  no  ripples.  However,  the  real 
physical  bottom  being  studied  should  exhibit  a  higher  roughness,  with  length 
scales  based  on  ripple  geometry  rather  than  the  sand  grain  diameter,  indi¬ 
cating  that  the  initially  chosen  value  of  =  0.01  is  probably  more  correct 
on  physical  grounds. 

In  this  regard,  we  note  the  questionable  practice  of  using  distinctly 
bi-directional  wave  data  to  generate  monochromat ic  input  wave  conditions  for 
calibration  purposes.  The  calibration  obtained  here  possible  constitutes  a 
valid  site  specific  calibration  for  the  Torrey  Pines  beach,  since  it  was 
found  to  be  possible  to  predict  net  longshore  flows  resulting  from  wave  fields 
with  several  components.  In  a  practical  sense,  few  data  sets  exist  which  are 
strictly  useful  for  calibration  purposes. 
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The  values  of  the  offshore  mixing  coefficient  N  chosen  in  this  study 
are  approximately  one  order  of  magnitude  smaller  than  the  corresponding  value 
suggested  by  Bowen  and  Inman  (1974),  who  investigated  a  group  of  dye  dispersion 
studies  performed  by  various  investigators.  Some  indication  that  a  larger 
value  of  N  may  be  desirable  is  given  by  Figure  (4-9):  however,  the  validity 
of  the  field  bathymetry  is  suspect  in  this  case.  It  should  be  noted  that 
the  model  exhibits  a  certain  degree  of  numerical  diffusion  when  large  grid 
spacings  are  used,  indicating  that  the  value  of  N  chosen  should  be  less  than 
the  corresponding  physically  realistic  value  in  any  case. 


Chapter  V 


EXAMPLES  OF  NUMERICAL  RESULTS 
USING  THE  NEARSHORE  CIRCULATION  MODELS 


1.  INTRODUCTION 

Various  results  of  calculations  using  the  linear  and  nonlinear  models 
have  been  described  in  detail  in  Birkemeier  and  Dalrymple  (1976)  and  Ebersole 
and  Dalrymple  (1979).  Results  pertaining  to  the  questions  of  model  stability 
and  convergence  have  been  mentioned  in  the  previous  chapters,  with  special 
attention  to  the  seiching  mode  of  oscillation  generated  at  model  start-up. 

In  this  chapter,  results  of  model  calculations  in  specific  situations  will 
be  discussed,  with  emphasis  on  comparison  to  analytic  models,  and  to  bottom 
topographies  which  reproduce  earlier  efforts. 

In  addition  to  the  general  facilities  for  user-defined  input,  each 
of  the  models  described  in  the  two  previous  studies  contained  specialized 
facilities  for  the  purpose  of  illustrating  specific  situations  not  covered 
by  the  basic  model  structures.  These  special  cases  are  discussed  here  for 
completeness,  although  in  most  cases  the  final  model  versions  may  not  retain 
the  corresponding  capability. 

2.  SPECIFIC  APPLICATIONS  OF  THE  LINEAR  AND  NONLINEAR  MODELS 

In  this  section,  we  review  results  presented  by  Birkemeier  and 
Dalrymple  (1976)  and  Ebersole  and  Dalrymple  (1979).  In  addition,  we  present 
results  for  an  additional  form  of  generalized  topography. 
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2 . 1  Plane  Beach  Applications 

Historically,  the  first  application  of  the  averaged  momentum  flux 
formulation  presented  in  Chapter  II  to  nearshore  dynamics  was  made  in  an 
attempt  to  explain  the  phenomenon  of  wave  set-up  and  longshore  wave-induced 
currents  in  the  surf  zone.  In  the  simplest  formulation,  the  equation  lead 
to  a  linear  longshore  current  profile  as  shown  in  Figure  2-2  (Longuet-Higgins, 
1970a).  The  addition  of  a  turbulent  stress  term  representing  the  effect  of 
lateral  mixing  leads  to  a  smoothed  profile  which  is  more  representative  of 
observed  current  distributions,  as  discussed  by  Longuet-Higgins  (1970b). 

In  the  development  in  Chapter  II,  a  lateral  mixing  model  in  the  x  and  y 
directions  was  outlined  based  on  the  model  used  by  Longuet-Higgins .  Here, 
we  review  results  calculated  both  with  and  without  the  lateral  mixing  effects. 

Birkemeier  and  Dalrymple  (1976)  presented  results  for  set-up  resulting 
from  normally  incident  waves  on  a  plane  slope,  as  shown  in  Figure  5-1  in 
comparison  with  the  latoratory  measurements  of  Bowen  et  al.  (1968).  Since, 
for  this  case,  no  currents  are  generated,  the  distinction  between  the  linear 
and  nonlinear  models  are  inapplicable.  The  model  is  seen  to  accurately  repro¬ 
duce  aspects  of  the  solution  based  on  linear  incident  waves,  including  the 
sharp  drop  in  the  mean  water  level  just  outside  of  the  breaker  line.  It  is 
noted  that  this  drop  can  be  eliminated,  and  a  more  realistic  solution  obtained, 
by  the  use  of  radiation  stress  terms  derived  from  cnoidal  wave  theory  (James 
1974)),  This  possible  refinement  has  not  been  included  in  the  models. 

Birkemeier  and  Dalrymple  also  investigated  the  dynamic  set-up  resulting 
from  a  time  varying  incident  wave  height  given  by 

H  =  H  +  A  sin  (27m  (5.1) 

os  1 

g 
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Figure  5-1.  Set-up  on  a  Plane  Beach 

12  -  Bowen  et.  al.  ( 1968)  experiment 
.  -  f, inear  Model 
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W 


where 


H  =  starting  wave  height 
A  =  amplitude  of  variation 
n  =  iteration  index 
T  =  group  period. 

Results  for  dynamic  set-up  in  comparison  to  incident  wave  height  are  shown 
in  Figure  5-2.  \s  an  example  of  the  seiching  response  of  the  model,  a  case 
was  also  run  where  the  group  period  corresponded  to  the  seiche  period  of  the 
model  as  given  by  Eq.  (3.5).  Figure  5-3  indicates  the  resulting  instability 
in  the  set-up  which  represents  the  growth  of  the  seiche. 

In  the  previous  example,  a  model  seiche  appeared  as  a  forced  response 
to  the  applied  wave-induced  stress.  Figure  5-4  indicates  a  more  typical  result 
of  a  seiching  response  to  the  transient  model  start-up.  In  this  case,  the 
oscillation  is  a  free  motion  which  gradually  dies  away  due  to  the  effect  of 
bottom  friction  (as  well  as  numerical  dispersion).  As  mentioned  in  section 
3.4,  this  type  of  response  can  be  reduced  by  a  gradual  build-up  of  the 
incident  wave  height. 

Results  for  a  longshore  current  on  a  plane  beach  are  shown  in 
Figure  5-5,  for  the  linear  model  and  no  lateral  mixing,  and  in  Figure  5-6, 
for  the  nonlinear  model  with  mixing.  The  two  solutions  clearly  represent 
the  behavior  predicted  by  the  theoretical  solutions,  with  the  exception  that 
the  sharp  decay  in  velocity  at  the  breaker  line  predicted  by  the  solution 
with  no  mixing  is  spread  over  several  grid  spaces  in  the  linear  model  duo 
to  the  use  of  finite  grid  squares.  This  effect  can  be  reduced  by  the  use 
of  a  finer  grid  spacing. 
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Figure  5-5.  Longshore  Current  Profile  for  Plane  Be 
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Figure  5-6.  Longshore  Current  Profile  for  Plane  Beach  including 

Lateral  Mixing 
Nonlinear  Model 


The  applications  discussed  here  indicate  that  the  model  accurately 
reproduces  the  available  analytic  solutions  for  simplified  physical  situations. 

2 . 2  Barred  Profile  Applications 

In  nature,  beaches  are  often  fronted  by  continuous  or  fragmented 
longshore  bars.  In  order  to  model  this  situation,  the  nonlinear  model  was 
run  using  a  barred  profile  (Figure  5-7)  for  the  same  wave  conditions  used 
to  generate  the  plane  beach  results  described  above.  For  this  case,  the 
model  was  run  neglecting  the  effect  of  lateral  mixing,  and  including  the 
effect;  results  are  shown  in  Figures  5-8  and  5-9  respectively.  The  solution 
without  lateral  mixing  demonstrates  the  major  discrepancy  found  between 
theoretical  solution  and  field  data  for  the  special  case.  As  the  model 
wave  passes  the  crest  of  the  offshore  bar,  it  stops  breaking  since  it 
responds  instan taneously  to  the  local  bottom.  In  the  absence  of  lateral 
mixing,  no  currents  are  generated  in  the  trough  between  shore  and  bar.  This 
result  is  in  disagreement  with  field  observation  as  shown  recently  by  Allender 
et  al.  (1981)  as  well  as  other  investigators.  The  inclusion  of  lateral  mixing 
effects  partially  alleviates  this  discrepancy. 

Two  factors  may  contribute  to  the  observations  that  currents  in 
natural  surf  zones  tend  to  be  strongest  in  the  trough  between  bar  and  shore. 
First,  fragmented  bars  tend  naturally  to  induce  two  dimensional  flow  patterns 
characterized  by  rip  currents  flowing  seaward  at  the  gaps  in  the  bar  system. 
These  currents  are  driven  by  longshore  variations  in  the  set-up  resulting 
from  interaction  between  the  incident  waves  and  the  non-uniform  topography. 
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BARRED  PROFILE 

PLANE  BEACH( SLOPE = .025 ) 


Figure  5-7.  Barred  Profile  used  for  Results  in  Figures  5-8,  5-9. 
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5-8.  Longshore  Current  for  Barred  Profile 
Nonlinear  Model,  No  Lateral  Mixing 


Figure  5-9.  Longshore  Current  for  Barred  Profile 
Nonlinear  Model,  Lateral  Mixing  Included 


89 


In  these  situations,  currents  flowing  along  the  beach  are  driven  by  gradients 


in  the  mean  water  surface;  the  flow  is  naturally  higher  in  the  trough,  where 
hydraulic  resistance  is  lower.  This  mechanism  for  rip  current  maintenance  has  been 
discussed  by  Dalrymple  (1978).  It  is  possible  that  this  mechanism  will 
always  bias  the  currents  observed  on  real  beaches. 

The  second  possibility  perhaps  lies  in  the  standard  treatment  of 
surf  zone  wave  height  as  a  function  only  of  the  local  depth.  The  process 
of  wave  breaking  initiates  a  turbulent  flow  pattern  which  is  certain  to 
have  some  time  dependence,  at  least  in  its  decay.  A  model  for  spilling 
breakers  including  this  effect  has  been  discussed  by  Longuet-Higgins  and 
Turner  (1974).  It  is  possible,  then,  that  an  accurate  picture  of  currents 
on  a  barred  profile  would  require  the  inclusion  of  a  continuation  of  wave 
energy  decay  in  a  time  dependent  manner  past  the  roint  where  the  wave  stops 
shoaling.  Wave  energy  decay  models  have  been  used  successfully  by  Divoky, 

LeMehaute  and  Lin  (1970)  in  a  study  of  wave  height  decay  in  the  surf  zone,  and 
by  Miller  and  Barcilon  (1978)  in  a  study  of  rip  currents.  However,  present 
models  do  not  allow  for  the  cessation  of  breaking,  which  is  required  if  the 
process  of  wave  breaking  and  reformation  is  to  be  successfully  treated. 

2 . 3  Applications  to  the  Laboratory  Wave  Basin 

In  order  to  model  currents  in  a  laboratory  wave  basin,  the  linear 
model  described  here  was  modified  to  include  no  flow  boundary  conditions  at 
the  longshore  boundaries.  (This  option  is  not  included  in  the  present  model.) 
Model  tests  were  conducted  to  numerically  approximate  the  experimental  set-up 
and  analytic  theory  of  Dalrymple  et  al.  (1977).  In  the  physical  experiment, 
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a  plane  beach  was  established  at  an  angle  of  15°  to  a  flap-type  wavemaker. 


Waves  approached  the  beach  over  a  flat  bottom  until  they  reached  the  front 
of  the  beach  slope,  where  refraction  effects  began.  Three  cases  were  tested 
experimentally  and  theoret ically ;  here,  we  restrict  our  attention  to  the  case 
where  the  surf  zone  is  bounded  laterally  by  walls  extending  to  infinity  in 
the  offshore  direction.  Experimentally  determined  streamlines  are  shown  in 
Figure  (5-10),  in  comparison  to  the  analytic  solution.  In  the  numerical 
model,  the  physical  situation  was  altered  by  extending  impermeable  walls 
in  a  direction  normal  to  the  beach;  the  waves,  however,  are  allowed  to  propa¬ 
gate  freely  according  to  refraction  governed  by  an  infinite  beach.  The 
numerical  model  corresponded  to  the  conditions  used  to  develop  the  analytic 
solution.  Velocity  vectors  calculated  by  the  linear  model  are  shown  in 
Figure  (5-11),  in  mirror  image  to  the  geometry  in  Figure  (5-10).  Longshore 
current  velocities  calculated  by  the  linear  model  are  compared  to  experi¬ 
mental  values  in  Figure  (5-12)  .  The  model  is  seen  to  underpredict  currents 
in  comparison  to  the  analytic  and  experimental  results. 


2.4  Periodic  Bottom  Topography 

In  a  study  of  rip  currents  caused  by  submerged  on-offshore  channels, 
Noda  (1970)  developed  an  equation  for  the  depth  which  produces  a  channel 
oriented  at  an  angle  to  the  beach.  The  present  models  were  tested  using 
Noda's  periodic  bottom,  given  by 


Experimental  Strealines  in  the  Wave  Basin  [From  Dalrvmple  et  al.  (1977)] 


VELOCITY  (cm  /sec) 


Figure  5-12,  Predicted  and  Measured  Longshore  Velocities 
for  the  Closed  Basin  Test  Case. 

Linear  Model. 
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0.  . 
l,J 


-m(4  -  i)Ax  i  =  1,2, 3, 4 

^mAx{  1+A  exp  [-  3(~)^^jsin^  y  (y-x  tan  B)}  i  >  4 


where  x,  v  are  coordinates  of  grid  (i,j) 

m  =  average  beach  slope  =  0.025 
A  =  length  of  periodicity  =  70  m 
A  =  amplitude  of  bottom  variation  =  20 
B  -  angle  of  rip  channel  to  beach  normal  =  30° 

The  test  contours  are  shown  in  Figure  5-13.  The  bottom  topography  was 
intended  to  model  a  field  case  studied  by  Sonu  (1972).  Current  vectors 
calculated  using  the  linear  model  are  shown  in  Figure  5-14;  current  magni¬ 
tude  agree  well  with  the  field  data  reported  in  Noda  £t  a_l,  (1974),  and  are 
higher  than  the  currents  calculated  in  Noda  et  al,fs.  model, 

in  which  only  50%  of  the  calculated  current  is  used  in  determining  the  wave 
number  due  to  numerical  instability  problems. 


The  Noda  profile  was  also  tested  using  the  nonlinear  model  with 
lateral  mixing  (Figure  5-15).  The  effect  o f  the  extra  terms  model  lei  is  to 
largely  drop  out  the  rip  current  clearly  seen  in  the  linear  model  results. 
Similar  results  have  been  obtained  by  Liu  (1982),  using  a  finite  element 
approach,  but  the  absence  of  a  rip  current  under  the  given  conditions  is 
clearly  at  odds  with  the  field  data  of  Sonu  (1972). 
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A  second  periodic  topography  has  been  investigated  for  the  present 
study.  It  is  well  established  that  the  presence  of  a  nearly  shore-normal 
channel  deeper  than  the  surrounding  beach  slope  will  tend  to  induce  a  rip 
current,  with  flow  directed  offshore  in  the  channel.  For  the  present  case,  a 
longshore-periodic  perturbation  consisting  of  a  localized  bulge,  or  relatively 
shallow  region,  has  been  added  to  an  otherwise  planar  beach  profile.  The 
bottom  is  given  by  the  relation 


D.  . 
I*  J 


f-  m  (2  -  i)  \x 


i  =  1 


(i‘x)(l  -  0.9  cos^  ^  (N^l)  ~0*015  i.'-x  )  i  1 


where 


m  -  average  beach  slope  =  0.02 

i,j  =  grid  coordinates 

(N-l).’.y  =  longshore  periodicity 


Results  using  both  models  are  shown  in  Figures  5-16  and  5-17.  The 
effect  of  the  inclusion  of  convective  acceleration  terms  in  the  nonlinear 
model  is  not  apparent  in  the  results,  which  may  be  a  result  of  the  counter¬ 
balancing  effect  of  lateral  mixing.  The  lateral  mixing  effect  can  be  seen 
in  the  shoremost  grid  row.  The  overall  effect  modelled  here  can  be  explained 
most  simply  in  terms  of  a  longshore  imbalance  in  the  steady  state  set-up. 

As  the  bulge  concentrates  wave  energy  by  refraction,  a  region  of  relatively 
high  breaking  waves  is  created.  As  water  is  pumped  on  shore  by  the  gradient 
of  the  onshore  radiation  stress,  a  localized  region  of  high  set-up  is  created; 
this  region  in  turn  pumps  the  onshore-flowing  water  alongshore  into  the 


regions  of  low  set-up.  The  waves  over  the  region  of  the  bulge  then  must 
continuously  pump  water  onshore  in  an  attempt  to  fill  the  surf  zone  up  to  the 


O'., 


36C 


320. OC 


280 . 


2C3.C 


160.00 


IOC. 00 


80  .OC 


0 


Figure  5-16. 


Currents  Induced  by  Bulge  on 
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a  Plane  Beach  Linear  Model: 
180°. 


100 

i 


Depth  below  Still  Water  Love]  (m) 


Depth  below  Still  Water  Level  (m) 


level  determined  by  the  decay  of  wave  energy,  leading  to  a  steady  state 
circulation  pattern. 

2 . 5  Intersecting  Waves  Application 

Ebersole  and  Dalrymple  have  discussed  the  application  of  the  model 
to  the  case  of  intersecting  wave  trains  of  the  same  frequency  on  a  plane 
beach,  which  provides  a  mechanism  for  generating  rip  currents,  as  shown  by 
Dalrymple  (1975).  The  purpose  of  this  application  was  to  investigate  the 
effect  of  the  convective  acceleration  terms  in  the  model.  The  following 
derivation  closely  follows  the  work  of  Dalrymple. 

Given  two  intersecting  wave  trains  A  and  B  with  amplitudes  a  and  b 
and  a  common  frequency,  a,  in  terms  of  the  coordinate  system  shown  in 
Figure  5-18,  the  free  surface  displacements  for  the  two  wave  trains  can 
be  written  as, 

nA  =  a  cos(k  cos  ax  4  k  sin  ay  4  at) 


n  =  b  cos(k  cos  Bx  4-  k  sin  By  4  at) 

D 


Figure  5-18.  Wave  Angles  for  Intersecting  Waves. 
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The  total  free  surface  n^,  =  can  then  be  written  as, 

h  k 

ri»j,  =  -a  cos  {~(cosa  +  cos  £)x  +  —  (sim  +  sin  3)y  +  ot} 
k  k 

cos  {— (cOvSa  -  cos  3)x  +  y(sina  -  sin  £)y) 


+  (b-a)cos  {k  cos  rx  +  k  sin  fv  +  at } 


(5.2) 


Using  the  linearized  dynamic  free  surface  boundary  condition  the  velocity 
potential  can  be  shown  to  equal, 

2ag  cosh  k (h+z)  .  fk/  ,  .  k,  .  _x  .  , 

^ — n — smivr(cosa  +  cos  3)x  +  7r(sina  +  sm  fc.)  +  it}  • 

T  o  cosh  kh  2  2 

c°s  (|(cosa  -  cos  ;)x  +  — (sinu  -  sin  S)y} 

,  (b-a)  cosh  k ( h+z )  r,  , 

+ - ; — n sln  Ik  cos  £x  +  k  sm  By  +  ct} 

a  cosh  kh  J 

From  the  velocity  potential  the  total  orbital  velocities  can  be  found  from, 


u  =  -  —  ,  v  — 

3x 


09 


H 


The  radiation  stresses,  which  are  essentially  the  forcing  terms,  are  defined 


as. 


r°  2 


xx 


pu  dz  + 


-h 


Pdz  -  j  pg(h-f-n)2  +  y  pg  n2 


yy 


xy 


0 

-h 

0 

-h 


pv  dz  + 


puv  dz 


-h 


Pdz  -  j  pg(h+n)2  +  -|  pg  n2 


where 


P  =  pg(n-z)  + 


3x 


puw  dz  -T 


Jy  J 


pvw  dz  -  pw" 
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The  radiation  stresses  are  found  to  be  given  by 


S 

xx 


_ tJL _ 

4sinh2kh 


2  2  2  2 
la  cos  a  +  b  cos  (3  +  2ab 


cosucosi:cos{  }  1  • 


{ 2kh+sinh2kh}  -  Q-P^-  ~~  (cosft-  cosa)4'cos{21^}  • 
8sinh2kh 

{  2khcosh2kh-sinh2kh  }  -  -  —  &  f  |?T  (sinx-sin;; )  "cos:  24> }  • 

8sinh2kh 


;  2khcosh2kh-sinh2kh I 


4sinh2kh 


[a‘'+b^+2ab  cos; 2ty}\ 


S  =  ; — oi  Z  [a^sin^a+b^sin^3+2absinnsin8cos{ 2^} ]  •  < 2kh+sinh2kh } 
yy  4sinh2kh 

-  qP~^  f v;  r  (cosp-cosa)^cos{ 2^}  •  { 2khcosh2kh~sinh2kh } 
8smh2kh 

-  f!?,  r  (sina-sinS) ^cosi 2^}  •  { 2khcosh2kh-sinh2khl 
8sinh2kh 

-  7 — «"P j'Tr r  [a^+b^+2abcos{ 2^}]  •  { sinh2kh-2kh } 

4sinh2kh 

2  1  2 
+  pgabcos  pg(b~a) 


S  =*  7 — .--f  n;  r  [a^sinacosa+b^cosBsin^+abcos{ 2 ^}sin  (a+p)  ]  * 

xy  4sxnhzkh 

{2kh+sinh2kh} 

where  the  expression  "ip"  is  defined  as, 

k  k 

^  =  -^(cosa-cosBjx  +  -j  (sina-sinp)  y 
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The  time  independent  mean  free  surface  displacement,  n,  is  defined 


where 


-  _  1  (  2^  2^  2  , 
n  -  “  -ZT-  iu  +V  +W  f 

2g  Z=0 


denotes  the  time  average  over  one  wave  period.  Substituting 


the  expressions  for  the  velocity  components  u,  v,  and  w  from  the  velocity 


potential  n  can  be  written  as, 

n  =  — : — nrrr  [a*"+b^+2abcosl2iJ;}#  (cos (a-£) cosh^kh-sinh^kh] 

2smh2kh 


(5.3) 


where  M^n  is  the  same  quantity  defined  previously.  Notice  that  the  mean 
free  surface  displacement  is  modulated  in  the  x  and  y  directions  by. 


cos{k(cosa-cos^)x  +  k (sinn-sinf )v } 


Using  Snellfs  Law  which  states 


k  sina  =  k  sin  ot  and  k  sinS  =  k  sin  j: 
o  o  o  o 


and  using  the  fact  that  k^  =  -j-“  *  where  f,on  denotes  deep  water  values  for 

o 

the  wave  length,  L,  and  the  wave  angles,  a  and  6,  we  see  that  there  is  a 
periodicity  of  the  mean  displacement  in  the  longshore  direction  with  a 


periodic  spacing,  given  by. 


sin  a  -  sin  3 
o  o 


This  periodicity  in  water  level  and  wave  height  causes  water  to  be  driven 
from  regions  of  high  mean  water  level  displacement  to  regions  of  lower  displace¬ 
ment,  resulting  in  the  formation  of  circulation  cells. 


103 


r 


In  order  to  attempt  to  model  tills  phenomena,  certain  simplifications 
to  the  model  had  to  be  made.  Since  the  refraction  and  shoaling  routines 
borrowed  from  the  work  of  Noda,ejt  al^.  ,  could  not  treat  more  than  one  wave 
train,  they  were  replaced  with  routines  governed  by  Snell’s  Law  neglecting 
wave-current  interaction.  Again  a  quadratic,  ’’exact,"  bottom  friction  was 
used  including  velocities  due  to  mean  currents  and,  this  time,  the  two  wave 
trains.  In  the  momentum  equations  the  advective  acceleration  terms  were 
retained,  horizontal  mixing  was  neglected,  and  the  radiation  stresses  were 
calculated  using  the  results  presented  earlier  in  this  section.  The  wave 
height  used  in  calculating  the  radiation  stresses  and  the  bottom  frictional 
stresses  is  given  by 

/?  2 

=  2.0  /a  +b  +2ab  cos  [k (cosa-cosu  )x  4-  k (sin ot-sin3)y ]  . 

Two  runs  are  presented  here  using  different  combinations  of  wave 
heights  and  wave  angles.  The  remainder  of  the  input  data  for  both  runs, 
however,  was  the  same.  The  waves  were  run  on  a  plane  beach  with  a  slope 
of  0.025.  The  planform  area  of  interest  was  comprised  of  25  grids  in  the 
x  direction  with  an  Ax  grid  size  of  5.0  meters,  and  21  grids  in  the  longshore 
direction  with  a  Ay  grid  size  of  4.0  meters.  The  time  step  was  chosen  to  be 
0.2  seconds  and  the  model  was  run  for  1500  iterations  for  both  cases.  A 
wave  period  of  7.159366  seconds  was  used,  which  resulted  in  theoretical  rip 
spacings  of  80.0  meters.  The  bottom  friction  factor  was  set  equal  to  0.12 
to  allow  the  system  to  reach  steady  state  after  the  1500  iterations  and  to 
decrease  the  magnitude  of  the  resultant  currents. 

The  first  case  used  waves  of  equal  heights  and  equal  angles  to  either 
side  of  the  beach  normal.  The  deep  water  wave  heights  were  0.25  meters  and 
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the  deep  water  angles  were  ±  30.0  degrees.  For  this  case,  referring  to 
Eq .  (5.3),  a- a  and  a=~B  resulting  in  a  free  surface  displacement  given  by 

ru  -  2a  cos(k  sin  a  y)cos(k  cos  *x  +  ">t) 

T  o  o 

This  free  surface  describes  a  wave  train  moving  in  the  -x  direction 
with  a  modulated  wave  height  that  is  periodic  in  the  longshore  direction 
only.  The  periodicity  in  wave  height  is  the  driving  mechanism  producing 
the  rip  current  perpendicular  to  the  beach,  as  shown  in  Figure  5-39.  Note 
the  constricted  width  of  the  rip  current  in  relation  to  the  width  of  the 
inflow  region.  This  is  a  result  of  the  convective  acceleration  terms. 

Also  note  the  weak  rip  head  where  the  currents  diverge  from  the  rip  axis 
and  return  towards  shore. 

In  the  second  case,  the  waves,  A  and  B,  had  different  heights,  0.1 
and  0.4  meters,  and  wave  angles  of  30.0  and  -30.0  degrees,  respectively. 

The  resulting  circulation  pattern  is  shown  in  Figure  5-20,  and  consists 
of  a  meandering  current  with  alternating  regions  of  strong  and  weak  long¬ 
shore  velocity  along  the  beach.  This  circulation  would  lend  itself  well 
to  the  formation  of  rythmic  beach  features.  Looking  at  Eq .  (5.2),  we  see 
that  there  is  a  non-zero  term, 

(b-a)cos(k  cos  (?x  +  k  sin  By  +  ct} 

which  is  a  wave  train  at  an  angle  to  the  beach  normal  with  height  2(b-a). 

This  wave  is  present  in  addition  to  the  normal  wave  train  with  the  modulated 
height  from  the  first  case,  causing  a  longshore  current  which  is  superimposed 
on  the  cellular  circulation. 
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Figure  5-19.  Current  Vector  Plot  for  Intersecting  Waves,  Case  1. 

Non li nea  r  Model 
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Figure  5-20,  Current  Vector  Plot  for  the  Meandering  Circulation  Pattern. 
Nonlinear  Model. 
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Chapter  VI 


CONCLUSIONS 

In  this  report,,  the  theoretical  background  and  numerical  formulation 
of  two  models  for  nearshore  circulation  have  been  reviewed.  A  review  of  the 
application  of  the  models  to  plane  beach  problems  has  shown  that  the  models 
successfully  reproduce  the  characteristics  of  known  analytic  solutions.  The 
comments  of  Basco  (1981),  who  expressed  the  opinion  that  the  models  are  in¬ 
valid  due  to  the  effect  of  numerical  viscosity  on  the  calculated  velocity 
profiles,  can  be  seen  to  be  incorrect  upon  inspection  of  Figure  5-5,  which 
demonstrates  the  steep  drop  in  longshore  current  at  the  breaker  line  in  the 
absence  of  a  prescribed  lateral  mixing.  In  addition.  Figure  5-1  shows  that 
the  model  reproduces  the  sharp  drop  in  wave  set  down  at  the  breaker  line, 
further  indicating  that  the  model  is  capable  of  reproducing  the  distinct 
features  of  analytic  solutions,  if  the  grid  scheme  chosen  is  fine  enough 
to  resolve  the  features.  It  is  apparent  that  any  numerical  diffusion  of 
solutions  for  n  or  V  is  confined  to  one  or  two  grid  rows,  and  can  be  made 
insignificant  by  choosing  the  grid  size  small  enough.  It  should  also  be 
remarked  that  the  slight  decrease  in  the  longshore  velocity  below  the  purely 
linear  profile  predicted  by  Longue t-Higgins  (1970a)  has  been  explaineu  to 
be  a  result  of  wave-current  interaction  (Dalrymple,  1980)  included  in  both 
models,  and  of  finite  angle  of  wave  incidence  (Liu  and  Dalrymple,  1978; 
Krauss  and  Sasaki , 1979) ,  included  in  the  nonlinear  model.  In  conclusion, 
there  appears  to  be  no  discrepancy  between  numerical  solutions  using  a  fine 
grid  and  the  corresponding  analytic  solutions. 
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In  practical  applications,  grid  spacings  are  often  chosen  which  lead 
to  some  false  numerical  averaging  of  the  solutions.  In  the  linear  model, 
this  averaging  manifests  itself  most  obviously  as  a  smoothing  of  the  long¬ 
shore  profile  across  the  breaker  line  in  the  same  manner  as  would  be  induced 
by  lateral  mixing  effects  (see  Figure  4-14).  This  indicates  that  the  mixing 
coefficients  included  in  the  nonlinear  model  should  be  smaller  than  expected 
on  physical  grounds  in  order  to  compensate  for  the  numerical  effects.  With 
the  choice  of  a  friction  factor  based  on  the  formulae  of  Jonsson  (1966)  or 
Kajiura  (1968),  the  model  is  then  seen  to  produce  realistic  current  profiles, 
both  in  form  and  magnitude,  for  the  field  site  chosen.  However,  it  should  be 
noted  that  no  detailed  data  set  currently  exists  which  satisfies  the  require¬ 
ment  of  a  monochromatic  input  x^ave  condition,  and  which  represents  a  field 
case  for  which  model  results  are  fairly  reliable,  such  as  the  absence  of  a 
large  longshore  bar. 

The  numerical  models  are  shown  to  be  successful  predictors  of 
nearshore  dynamics  in  situations  dominated  by  refraction  and  weak  wave- 
current  interaction  effects.  In  this  regard,  it  is  noted  that  the  model 
as  currently  formulated  cannot  handle  certain  effects,  such  as  the  influence 
of  finite  barriers  in  the  wave  field,  or  strong  interaction  of  waves  and 
opposing  currents.  The  modelling  of  both  of  these  effects  requires  the 
inclusion  of  a  capability  for  handling  wave  diffraction. 

Finally,  the  present  models  require  an  input  wave  condition  based 
on  a  monochromatic  wave,  whereas  wave  trains  in  nature  tend  to  include  a 
spectrum  of  waves.  In  the  case  where  the  incident  wave  spectrum  is  suffi¬ 
ciently  narrow-banded  in  frequency,  the  incident  wave  can  be  represented  as 
a  wave  of  single  frequency  with  a  modulated  amplitude;  such  an  effect  can 
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be  modelled  after  modification  of  the  programs  to  handle  a  time  varying 
deepwater  wave  height.  However,  modelling  of  broad-banded  spectra  at  the 
level  of  the  treatment  here  remains  a  complex  problem. 
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Appendix  I. 


USING  THE  NEARSHORE  CIRCULATION  PROGRAM 

The  nearshore  circulation  model  as  described  above  is  supplied  in 
two  versions,  referred  to  as  the  linear  and  the  nonlinear  versions.  Complete 
listings  of  the  programs  for  each  version  follow.  The  programs  as  written  are 
designed  to  run  on  a  Burroughs  B7700.  The  line 

$ RESET  FREE 

which  appears  on  the  front  of  each  program  allows  for  the  use  of  standard 
Fortran  on  the  B7700.  This  line  should  be  removed  before  operating  the 
program  on  a  different  system. 

The  two  versions  of  the  model  are  currently  designed  to  be 
indistinguishable  to  the  operator,  with  the  input  file  containing  wave 
and  wind  information  and  the  depth  grid  being  of  identical  format  for  each 
program.  The  output  file  generated  currently  at  a  line  printer  is  also 
identical  for  each  program.  The  exception  to  the  uniformity  between  the 
two  models  occurs  in  the  output  file  stored  at  the  end  of  execution  for 
subsequent  restart  of  the  model.  Since  the  nonlinear  version  of  the  program 
requires  the  storage  of  three  time  levels  of  data,  two  extra  time  levels  are 
included  in  the  output.  Thus,  the  nonlinear  version  of  the  program  can  not 
be  started  based  on  intermediate  results  generated  by  the  linear  version, 
as  not  enough  information  is  present.  It  is  possible  in  principle  to  start 
the  linear  version  of  the  program  using  intermediate  results  from  the  non¬ 
linear  version,  although  at  present  the  READ  statement  in  the  linear  version 
is  not  structured  to  handle  this  option. 
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The  circulation  model  is  currently  designed  to  be  run  as  a  batch 
job,  and  as  such  has  all  device  number  specifications  for  data  files  in  the 
job  file  external  to  the  program.  This  feature  may  have  to  be  changed  depending 
on  the  machine  to  be  used.  Currently,  data  files  may  be  named  arbitrarily.  The 
program  currently  requirues  four  10  device  specifications: 

Logical  Device  Number  Corresponding  Data  File 

5  User  defined  input  data 

6  Output  stream,  currently 

directed  to  a  line  printer 

3  Output  file,  stores  results 

of  last  interation  for  a 
later  restart 

8  Input  file  used  to  restart 

the  model 


It  should  be  noted  that  the  input  file  8  used  for  restarting  the  program 
is  identical  to,  and  should  be  just  a  renamed  version  of,  the  output  file  3 
generated  by  the  model  at  the  point  from  which  it  is  desired  to  continue 
computation. 


Structure  of  Input  Data  File  5 

The  input  data  file  is  structured  into  five  groups  of  data.  The  first 
four  groups  each  consist  of  a  single  line  in  the  file,  with  data  values  entered 
unformatted;  i.e.,  separated  by  commas  with  no  requirements  on  spacing.  The 
fifth  data  group  consists  of  the  depth  grid,  and  requires  more  space  than  a 
single  line,  as  explained  below.  The  information  included  in  the  input  data 
file  is  as  follows. 
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DATA  IN  INPUT  DATA  FILE  5 


I .  Wave  Pa rameters 

T  -  Wave  period  (seconds) 

HO  -  Wave  height  (meters) 

A  -  Wave  angle,  measured  clockwise  from  the  -fX  direction 

(degrees) 


2.  Wind  Parameters 

WIND  -  Wind  speed  (meters/sec.ond) 

WINANG  -  Wind  angle,  measured  clockwise  from  the  -X  direction 
(degrees) 


3.  Grid  Parameters 

M  -  Number  of  grid  rows  in  X  direction  (offshore) 

N  -  Number  of  grid  columns  in  Y  direction  (longshore) 

DX  -  Grid  spacing  in  X  direction  (meters) 

DY  -  Grid  spacing  in  Y  direction  (meters) 

DT  -  Size  of  the  time  step  (seconds) 

INDEX  -  Specifies  the  input  of  depth  information 

=1,  read  data  from  previous  run  from  file  8 

-2,  read  depth  grid  from  input  file  5 

=3,  establish  plane  beach  based  on  input  beach  slope 

AM  -  Beach  slope  (space  in  input  file  must  be  filled,  but 

value  is  unsed  only  if  INDEX=3) 


4.  Program  Control  Parameters 

ITA  -  Total  number  of  iterations  (including  the  accumulated 

total  from  previous  runs,  if  previous  run  data  is  used) 

NHIGHT  -  Number  of  iterations  over  which  the  deep  water  wave 
builds  up 


121 


ID  -  Determine  if  dissipation  of  wave  energy  is  to  be 

included 

=0,  no  dissipation 

-1,  dissipation  due  to  viscosity  and  bottom  permeability 
is  included 

KSKIP  -  Determine  frequency  of  printed  iterations  in  output 

file  6 

5.  Depth  Grid  (required  if  INDEX=2) 

D  -  Local  water  depth  with  respect  to  still  water  level 

at  each  grid  center. 

The  depth  grid  D  is  input  in  an  unformatted  string  of  length  M*N, 
starting  with  the  fifth  Line  of  the  input  file  5  and  continuing  to  as  many 
lines  as  necessary.  The  data  values  are  read  off  of  an  established  grid 
row-wise  in  the  +Y  direction,  starting  at  the  inshore  end  of  the  grid. 

Coordinate  System  Used  in  the  Program 

The  input  depth  grid  and  wave  and  wind  parameters  should  be  established 
using  the  coordinates  shown  in  Figure  2-1.  The  models  require  that  the  bathy¬ 
metry  be  periodic  in  the  longshore  direction.  Accordingly,  for  an  input  depth 
grid  of  size  (M,  N) ,  where  (I,  J)  are  X  and  Y  indices  and  (M,  N)  are  the  maximum 
values  of  (I,  J),  the  laterally  bordering  depth  values  should  satisfy  the 
requirement : 


D(I,N)  =  D(I,1) 

If  this  condition  is  not  met  on  input,  D(I,1)  is  redefined  to  be  equal  to  D(I,N). 
Both  models  create  two  additional  rows  in  the  Y  direction,  (N+l)  and  (N+2) .  The 
grid  storage  requirement  for  all  variables  is  thus  M  by  (N+2),  with  a  current 
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maximum  of  50  by  50.  The  maximum  may  be  enlarged  or  decreased  by  alteration 
of  all  COMMON  and  DIMENSION  statements  in  either  program. 

Wave  angle  and  wind  angle  are  defined  as  shown  in  Figure  2-1.  Wave 
angle  is  measured  clockwise  from  the  +X  direction.  Waves  approaching  the 
beach  from  the  right  of  directly  offshore  have  wave  angles  >  180°.  Wind 
angle  is  given  clockwise  from  the  -X  direction,  so  that  winds  approaching 
shore  from  the  right  of  normal  have  wind  directions  greater  than  0°. 

The  program  requires  that  incident  waves  have  longshore  components 
propagating  in  the  +Y  direction.  For  the  case  of  wave  angles  =  180°  -  6,  the 
program  flips  the  depth  grid  over  and  rotates  the  wave  and  wind  directions 
through  360°  -  20  degrees,  and  runs  the  program  for  waves  approaching  at 
180°  +  0  and  a  mirror  image  depth  grid. 

Error  Messages  from  the  Circulation  Models 

Each  version  of  the  nearshore  circulation  model  will  generate  various 
error  messages  when  fatal  conditions  are  met  during  execution.  The  messages 
are,  for  the  most  part,  identical  from  either  program  and  correspond  to  equivalent 
situations  within  the  programs.  A  list  of  possible  internally  occurring  errors 
leading  to  the  printing  of  a  message  follows. 

1.  Failure  of  iteration  scheme  to  converge  on  a  wave  height  at  a 
grid  point  will  cause  the  message 

RELAXATION  FOR  THE  WAVE  HEIGHT  FAILED  TO  CONVERGE 
followed  by  an  indication  of  the  row  and  number  of  iterations  tried.  This 
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condition  arises  in  the  subroutine  HEIGHT.  Occurrence  of  the  condition  does 
not  lead  to  termination  of  the  program;  iteration  continues  with  the  last  wave 
height  calculated  and  the  entire  iteration  is  tried  again  at  the  next  time  step. 

2.  Failure  of  iteration  scheme  to  converge  on  a  wave  angle  at  a 
grid  point  will  cause  the  message 

RELAXATION  FOR  THETA  FAILED  AFTER  —  ITERATIONS 

with  an  indication  given  of  the  number  of  iterations  tried.  This  condition 
could  arise  if  the  magnitude  of  a  current  directed  against  the  general  direction 
of  wave  travel  was  too  strong.  The  condition  arises  in  the  subroutine  ANGLE 
and  causes  termination  of  the  program. 

3.  If  a  negative  wave  number  RK  is  calculated,  the  message 

RK  IS  NEGATIVE 

is  printed,  and  the  program  is  terminated.  This  condition  occurs  in  the 
subroutine  WVNUM. 
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